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|  INTERACTION  BETWEEN  WATER  BORNE  WAVES  ANO  SEISMIC  WAVES  IN  THE  OCEAN 

[>  BOTTOM:  THE  FORWARD-  ANO  INVERSE  PROBLEM. 

ABSTRACT 

I  The  work  is  divided  into  two  mam  parts.  The  first  part,  the 

forward  problem,  is  aimed  at  establishing  a  theoretical 
,  framework  for  excitation  and  propagation  of  elastic  waves  in 

linear  homogeneous  isotropic  media. 

The  second  part,  the  inverse  problem  ,  is  aimed  at  determining 
the  environmental  parameters  which  significantly  influence 
acoustic  propagation  in  a  shallow  water  environment. 

The  methods  developed  here  indicate  that,  under  favorable 
conditions,  it  is  possible  to  infere  ocean  bottom  parameters 
sucn  as  P-  and  S-  wave  phase  velocities  from  near  surface 
measurements . 


I  INTRODUCTION 

The  work  reported  here  is  a  partial  fulfillment  of  the  requirements 
for  the  degree  Dr  Ing  under  the  supervision  of  professor  Jens  Hovem  at 
the  University  of  Trondheim,  the  Norwegian  Institute  of  Technology. 

The  work  is  divided  into  two  mam  parts.  The  first  part,  the  forward 
problem,  is  aimed  at  establishing  a  theoretical  framework  for 
excitation  and  propagation  of  elastic  waves  in  linear  homogeneous 
isotropic  media.  As  the  solution  is  limited  to  horizontally  stratified 
media,  a  separable  solution  of  the  equation  of  motion  is  derived. 

In  the  outset,  a  numerical  solution  was  implemented  based  on  the 
Thomson-Haskel  matrix  method[T,4]  as  derived  by  Kutschalel  T]  .  This 
approach  was  late-  abandonded  and  a  two-dimensional  model  -  developed 
at  the  SACLANT  ASW  Research  centre  at  LaSpezia,  Italy  by  Henrik 
Schmidt  [9]  -  was  modified  to  the  present  three-dimensional  version. 


As 
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The  second  part,  the  inverse  problem  ,  is  aimed  at  determining  the 
environmental  parameters  that  significantly  influence  acoustic 
propagation  in  a  shallow  water  environment.  This  is  done  by 
interpreting  shot  data  gathered  by  a  seismic  survey.  Interpretation  is 
sought  verified  by  comparison  with  output  from  the  numerical  model.  To 
my  knowledge,  the  inverse  problem  methods  utilized  here  have  not  been 
published  earlier. 
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7  THE  FORWARD  PR08LEM 

2 . t  Introduction 

It  is  well  known  that  separable  solutions  aad  to  the  use  of  integral 
transform  techniques  which  yield  an  exact  solution  to  the  wave 
equation  in  stratified  elastic  mediatl],  The  field  parameters  are. 
however,  determined  by  linear  combinations  of  the  basis  for  the 
solution  space  i  e  inverse  transform  integrals.  In  cases  with  only  a 
few  layers,  contour  integration  can  be  used  to  reduce  the  numerical 
computation.  involving  only  a  few  integrations  over  finite  intervals, 
e  g  [1],  In  general  numerical  models.  however,  such  techniques  are 
inconvenient,  and  direct  numerical  integration  has  to  be  used. 

In  underwater  acoustics  in  general,  and  also  for  the  cases  we  will 
investigate  here,  the  sources  are  usually  contained  within  a  volume 
small  compared  to  the  volume  of  interest,  thus  the  radiated  field  is 
most  conveniently  described  in  a  cylindrical  coordinate  system.  The 
field  is  then  given  by  Hankel  transform  integrals  which  are  not  well 
suited  for  direct  numerical  integration  due  to  the  Bessel  functions 
involved.  In  order  to  overcome  this  problem,  Marsh[2]  in  19B1 
introduced  what  was  later  called  the  fast  field  approximation  IFFP  = 
Fast  Field  Programl  of  the  Hankel  transform.  The  field  is  separated 
into  ingoing  and  outgoing  parts  by  expressing  the  Bessel  function  in 
terms  of  Hankel  functions,  the  ingoing  part  is  disregarded  and  the 
outgoing  part  is  replaced  by  its  large  argument  approximations.  The 
integrals  are  then  evaluated  by  means  of  the  fast  Fourier  transform 


( FFT  )  .  As  shown 

later 

by 

DiNapoli 

and 

Deavenport[5] , 

for  the 

two-dimensional 

case , 

the 

fast 

field 

approximation 

gives  no 

significant  errors  at  ranges  longer  than  a  few  wavelengths  from  the 
axis  . 

After  the  introduction  of  the  fast  field  technique,  a  number  of 
numerical  models  have  been  developed,  based  on  this  integration  method 
and  thus  usually  called  fast  field  programs. 
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In  spite  of  their  common  name,  these  models  are  significantly 
different,  especially  concerning  the  approach  taken  to  solve  the 
transformed  wave  equations  in  a  multilayered  environment. 
Traditionally,  the  depth  dependence  of  the  field  has  been  determined 
by  means  of  the  Thomson-Ha skell  matrix  method[3,4]  .  The  first  model 
was  introduced  by  DiNapoli[6],  who  evaluated  the  solution  very 
efficiently  by  means  of  recurrence  relations  for  the  hypergeometric 
functions.  However  this  approach  allows  only  for  fluid  layers,  and  in 
that  case,  other  techniques,  like  normal  mode  methods,  are  usually 
more  convenient.  The  first  FFP  model,  including  the  coupling  between 
compres s lonal  and  vertical  shear  waves  at  the  boundaries  of  solid 
layers,  was  developed  by  Kutschale[7]  also  using  the  Thomson-Haskell 


method.  The 

original 

model 

allowed 

for 

only 

one 

source/ receiver 

combination 

for 

each 

solution.  It 

has 

later 

been 

modified 

by 

Harrison [ 8 ] 

to 

allow 

for 

several 

receivers , 

but 

even  for 

one 

combination  the  computations  are  rather  extensive. 

A  more  direct  and  computationally  more  efficient  solution  technique 
was  recently  introduced  by  Schmidt[9].  The  field  parameters  at  the 
interfaces  are  expressed  in  terms  of  source  contributions  and  unknown 
scalar  potentials.  The  boundary  conditions  yield  a  system  of  equations 
in  the  Hankel  transforms  of  the  potentials  to  be  satisfied  at  each 
interface.  These  local  systems  of  equations  are  mapped  into  a  global 
set  of  equations  using  a  technique  similar  to  the  one  used  in  finite 
element  programs.  The  computational  speed  has  been  improved  by  an 
order  of  magnitude  by  use  of  this  solution  technique.  Furthermore 
configurations  involving  several  sources  and/or  receivers  can  be 
treated  with  one  solution,  thus  yielding  the  possibility  of  computing 
total  fields  generated  not  only  by  single  point  sources,  but  also  by 
vertical  source  arrays. 

These  and  similar  models  have  all  been  two-dimensional,  thus 
restricting  the  sources  to  be  placed  on  the  axis  of  the  cylindrical 
coordinate  system.  A  direct  solution  of  problems  with  horizontally 
distributed  sources  has  therefore  not  been  possible,  but  has  required 
a  new  calculation  for  each  source  and  subsequent  superposition.  In 
this  paper  the  model  of  Schmidt[9]  has  been  modified  and  extended  to 
allow  for  sources  displaced  with  respect  to  the  axis. 
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The  field  parameters  are  expanded  in  a  Fourier  series  in  the  angular 
direction,  thus  leading  to  an  infinite  number  of  two-dimensional 
problems.  By  expressing  the  boundary  conditions  in  terms  of  Cartesian 
components,  rather  than  polar  components,  the  coefficient  matrix  will 
be  independent  of  the  Fourier  order,  and  the  Hankel  transforms  of  all 
the  expansion  coefficients  for  the  unknown  potentials  can  be  found 
with  only  one  matrix  inversion  for  each  horizontal  wavenumber.  The 
truncation  point  of  the  Fourier  series  can  be  determined  a  priori. 

The  inversion  of  the  Hankel  transforms  is  again  performed  by  means  of 
the  fast  field  technique,  and  the  angular  distribution  is  evaluated 
from  the  expansion  coefficients  by  means  of  an  FFT  technique. 

In  the  following  the  model  and  its  mathematical  background  will  be 


described . 
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2.2  List  of  symbols 


hi 

A 

M 

Q 

c 

c 

c 

s 

<p 

* 

A 

-* 

0 

u 

u 

l 

-♦ 

V 

v 

1 

T 

Ik 


S 

c 


lk 

lklm 


tk 

5(kI 


V*’ 

N  Iz) 
h|Z1 (Z) 


K 

0. 


angular  frequency 

Lame'  constanc  and  wavelength 

tame'  constant 

density 

congressional  wave  phase  velocity 

shear  wave  phase  velocity 

compressional  wave  potential 

shear  wave  potential 

shear  wave  potential 

shear  wave  potential  vector 

particle  displacement  vector 

particle  displacement  in  l-direction 

particle  displacement  velocity  vector 

particle  displacement  velocity  in  i-direction 

stress  in  direction  i  applied  to  surface  with  normal 

vector  k 

strain  in  direction  l  of  surface  with  normal  vector  k 

stiffness  constants 

Kronecker  delta 

Oirac  delta  "function" 

Bessel  function  of  first  kind,  order  l 
Neumann  function  of  order  1 
Hankel  function  of  second  kind,  order  1 
horizontal  wave  number 

vertical  wave  number  for  compressional  waves 
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p  vertical  wavenumber  for  shear  waves 

tilde:  denotes  source  terms 

Fm  vector  of  parameters  entering  boundary  conditions  in 

n 

layer  n  for  angular  order  m 

B™  vector  of  up-  and  downward-  going  potentials  in  layer 

n  for  angular  order  m 

A  matrix  relating  B™  to  Fm  m  layer  n,  i=l:  lower 

n .  i  n  n 

i=u:  upper 

m  ... 

R  source  contributions  in  layer  n,  agular  order  m 

n 
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2.3  The  homogeneous  solution 

The  solution  will  be  restricted  to  a  homogeneous  medium,  hence  the 
environmental  parameters  are  independent  of  spatial  position. 

Newtons  second  law  of  motion  for  a  solid  may  be  stated  in  a  frame  of 
Cartesian  coordinates  using  tensor  notation  as 


3T  .  { x  , x  , x  ) 
ik  x  v  z 

dx. 


3  U  ( x  ,  x  ,  x  ) 
I  X  v  z 

at2 


<  1 1 


The  subscripts  take  on  values  from  1  to  3  indicating  directions  x.y 
and  z  respectively.  Einstein's  summation  convention  l  e  summation 
over  identical  indexes  is  assumed. 


Ti(<  is  the  stress  tensor:  stress  in  direction  i  applied  to  surface 
with  normal-vector  k. 

is  particle  displacement  in  i-direction. 

The  spatial  arguments  of  T  and  u  will  in  the  following  be  omitted,  but 
assumed  tacitly. 

p  is  the  density  of  the  medium. 

It)  thus  equates  the  force  in  direction  i  applied  to  an  inf mitisimal 
volume  element  to  its  inertia  in  direction  i. 

As  the  only  forces  acting  on  the  volume  element  are  the  internal 
volume  forces,  the  volume  is  considered  source-free. 

The  medium  described  by  (l)  is  furthermore  lossless.  Viscous  losses 
can  be  included  by  entering  a  first  order  time  derivative  term.  A 
discussion  on  the  causality  of  t imedependen t  solutions  when  entering 
losses  is  given  in  [IT], 
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The  linear  elastic  material  properties  for  a  general  solid  may  be 
stated  as 


ik 


c  ,  ,  S, 
lklm  lm 


(2) 


which  is  Hooke's  law.  The  deformation  is 
3u 


s  , 

lm  2  3x 

m 


3x, 


] 


13) 


We  have  the  following  symmetry  relations 


lk 

^lm 


"lklm 


ki 


ml 

c 


kilm 


'ikml 


"lmik 


IU) 

Ub) 

Uc) 


(la)  follows  as  a  consequence  of  the  torsional  moment  vanishing  at 
equilibrium  [16],  (lb)  is  an  obvious  consequence  of  (3).  The  three 
first  terms  of  Uc)  are  consequences  of  (la)  and  (lb)  respectively, 
while  the  fourth  term  follows  from  thermodynamical  considerations 
[16], 


Thus.  has  at  most  21  independent  components. 


For  c  ,,  to  be  isotropic,  it  is  required  that  [11] 
lklm 


C .  .  ,  -  A6  ,  5,  *  u8  ,  6,  ♦  v5  5,  , 

lklm  lk  lm  ll  km  lm  kl 


15) 


(la) impl i es  that  y  =  v  since  c  ,  ,  ■=  c,  , 

lklm  kilm 

The  constitutive  relation  for  an  isotropic  solid  is  thus  simplified  to 


be 

3u 

3u , 

3uk 

Tik 

-  6 .  A  7  ”  ♦ 

lk  6x 

+  3\ J 

i 

(6) 

which  rewritten  in  V-operator  notation  is 


T  ,  =  8  A  9  •  U  ♦  2  P  S  , 

lk  ik  ik 


(7) 
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The  first  two  symmetry  properties  given  in  It)  state  that  both  T  and 
S  are  symmetric.  If  they  furthermore  are  real,  they  may  be  written 
as  real  symmetric  matrixes  which  consequently  are  selfadgoint. 
Oiagonalization  gives  the  stress  and  strain  respectively  along  what  is 
termed  the  principal  axes.  Then  t!5] 


X  7.U  ♦  2  p  S 


This  form  is  valuable  for  giving  insight  into  the  physical 

interpretation  of  the  Lame'  constants.  For  a  fluid  where  p  =  0.  X 

becomes  the  bulk  modulus.  For  a  solid,  the  effect  due  to  the  principal 

stress  T  would  be  a  dillitation  pluss  an  extension  in  the  n-direction 
n 

It  should  be  noted  that  the  material  constants  which  normally  enter 
into  our  equations  are  the  adiabatic  terms  as  we  assume  that 
negligible  heat  is  transfered  to  the  medium. 


In  order  to  combine  equations  (1)  and  (6),  (6)  is  differentiated  with 
respect  to  x^  and  the  common  term  eliminated.  The  result,  rewritten  in 
7-operator  notation,  becomes 


(X*p)77-U  *  p7  u 


Use  of  the  vector  identity 
-»  -*  ?■* 

7xVxA  -  77"A  -  7  A 


enables  us  to  bring  (9)  to  a  commonly  used  form  termed  the  equation  of 
motion  for  a  solid  isotropic  homogenous  lossless  medium  containing  no 


c  77.  U 
c 


where  the  phase  velocities  for  compressional-  and  shear  waves  are 


112a) 
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S  =  e  l,2b) 

respectively . 

The  introduction  of  the  wave  velocities  at  this  stage  is  an 
anticipation  of  things  to  come.  It  will  be  shown  that  these  quantities 
are  phase  velocities  and,  indeed,  phase  velocities  related  to 
compressional-  and  shear-waves  respectively. 

Separating  the  solution  of  (11)  into  a  space-  and  a  time-dependent 
part,  will  imply  a  timedependent  part  of  the  form  exp(^gwt)  of  which 

exp(guit)  is  chosen.  The  symbol  U,  up  to  now  used  to  denote  the  time- 
and  space-dependent  particle  movement,  will  in  the  following,  unless 

otherwise  stated,  be  taken  to  mean  gust  the  space-dependent  part.  This 
should  not  cause  confusion.  (11)  is  thus  reduced  to 

2  -»  2  ■»  2  -»  .» 

Cc  9V<U  -  c^  7x7xU  ♦  w  U  =  0  113) 

A  common  approach  for  further  reduction  of  (13)  1 1 2  3  C 1 3  ]  is  to  express 
the  solution  as  the  sum  of  a  longitudinal-  and  a  transversal  term  i  e 


U  =  U, 


where 

v-uT  =  0 

VxUL  =  0 

and  (15b)  into  (13)  yields 
=  0 

We  define 


Insertion  of  (11),  (  15a ) 

2  ■*  2  •*  2 ■* 
C  V7*U,  -  c  VxVxU.  *  m  U 
C  L  s  T 


(U) 


(15a)  - 


(15b) 


(16) 
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*  • 

O 

(17a) 

-♦ 

-* 

Q 

-  YxU T 

(17b) 

which  inserted  into  (16)  yields 

2  2 

c  c 

U  =  -4  VxQ  M8) 

2  2 

ID  UJ 

The  vector  quantity  in  the  second  term  of  (18)  may  be  reduced  to  a 

scalar  by  expressing  the  potential  vector  as 
2 

—  ‘‘q  r  Ya  ♦  VxAa  (19) 

c 

s 

where  a  is  a  unit  vector.  It  is  common  practice  in  z-dependent  media, 
C31C4JC7]  to  choose  a  to  be  aligned  with  the  z-axis.  The  displacement 
components  arising  from  Y  will  not  have  components  in  the  z-direction, 
thus  representing  a  transversal  shear  wave  which  does  not  enter  into 
the  two-dimensional  (range  and  depth)  case. 

The  displacement  may  now  be  written 


U  =  *  YxYz  ♦  YxVxAz  (20) 

As  (16)'  is  a  linear  operator,  here  termed  L,  then 

LIU  I  »  L ( V*  *  VxYz  ♦  VxVxAz)  =  0  (21) 

is  equivalent  with 

L(7Y)  ‘  L(VxYz)  ♦  L (VxVxAz )  =  0  (22) 

As  the  left  hand  side  of  (20)  is  obviously  a  vector  function  of  the 
three  spatial  directions,  it  may  be  decomposed  into  three  linearly 
independent  scalar  functions.  So  must  also  than  be  the  case  with  the 
right  hand  side  of  (20!  as  well,  the  three  scalar  functions  *,Y  and  A 
are  linearly  independent.  As  each  term  of  (22)  is  linearly 
independent,  they  can  at  most  be  a  constant.  We  will  here  in  the 
following  restrict  our  attention  to  assuming  this  constant  to  be  zero. 
Thus.  (22)  is  implied  by 
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UV*1  =  0  123a) 

L(9xYz)  =  0  123b) 

L ( VxVxAz )  =  0  123c) 

thus  (16)  is  satisfied  if 

Vlc2V2*  *  w2*)  -  0  124a) 

c 

Vxl-c^xVxY  .  u2'V )  =  0  (24b) 

s 

2  2 

Vx7x(-C  7x7xA  ♦  u>  A 1  =  0  (24c) 

s 


The  quantities  within  brackets  of  (24a)  can  at  most  be  constant.  As 
sources  at  infinity  are  discarded,  the  radiation  condition  imposes 
that  ♦  -*0  when  Irl-*",  thus  implying  the  bracketed  term  been  g  0.  (24b) 

and  (24c)  are  clearly  fullfilled  if  the  bracketed  terms  are  constant. 
The  same  physical  considerations  as  applied  to  (24a)  are  valid,  thus 
yielding  that  the  bracketed  terms  vanish  when  |’r|-»-. 

Application  of  the  vector  identity  (10)  to  the  bracketed  terms  of 
(24b)  and  (24c)  bring  these  two  equations  to  beeing  of  the  Helmholz 
type  and  in  summarizing:  if  (20)  is  a  solution  of  (13),  then 


V‘  * 

*  h2 

*  =  0 

V2  T 

*  g2 

r  =  o 

7 2  A 

2 

*  g 

A  =  0 

wher  e 

h  and 

g  are  the 

wavenumber  3 

respectively  : 

2 

2 

O) 

Q 

h  = 

2 

A 

♦  2u 

C 


(25a) 

(25b) 

(25c! 

for  compr es s lona 1  -  and  shear  waves 

(26) 


U 


(27) 
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It  has  been  stated  (11).  without  beeing  shown,  that  it  is  the 
compressional  and  shear  wave  velocities  that  enter  into  (9).  A 
separable  solution  of  (25)  in  a  Cartesian  coordinate  system  is 


0(x,y,z) 


-](kx«ky*kz! 
e  x  y  z 


where 


h 


2 


z 


(28) 


(29) 


As  all  directions  are  equivalent,  propagation  in  the  x-direction  may 
be  considered  with  no  loss  of  generality.  Inclusion  of  the 
timedependent  factor  yields 


♦  1 1 ,  x ) 


i(uit  -  k  x) 
e  x 


(30) 


For  constant  phase  the  argument  of  the  exponential  function  is 
constant,  thus 


x  ui  ui  A  ♦  2u  .  1  /  2 

-  =  —  =  -=(  )  =  c 

t  k  h  g  c 

x 


(31  ) 


Similarily,  a  solution  of  (25b)  and  (25c)  will  yield 


■M-  .  (  U  ,  1/2 


(  32  ) 


Application  of  (20)  to  (28)  gives 


-]kx 

j(uit  - 
e 

k  x ) 

X 

-• 

X 

(33a) 

■]kx 

j  C  uit  - 
e 

k  x) 

X 

-• 

y 

(33b) 

k2 

X 

j  (u»t  - 
e 

X 

X 

X 

>* 

2 

(33c) 
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♦  represents  a  wave  propagating  in  x-direction  with  phase  velocity  c^ 
and  with  particle  displacement  in  the  direction  of  propagation,  while 
¥  and  A  propagate  in  the  same  direction  with  phase  velocity  c^  but 
with  particle  displacements  perpendicular  to  the  propagation 
direction . 


If  U  =  (u  ,u„  ,u  }  are  the  displacement  components  referred  to  a 
n  rn  Bn  zn 

cylindrical  coordinate  system  {r,9,zi  in  a  homogeneous  layer  n.  113) 
will  be  satisfied  for  (20)  expressed  in  cylindrical  coordinates  as 


3*  1  3¥  32A 

_ n  _ n  _ n 

urn"  3r  r  38  *  3r3z 


(31) 


i  3* 

dY 

i  32a 

135) 

_  _ n 

_ n 

n 

r  38 

dr 

r  383z 

3*  13  3  1  32 

7*  -  1  -  ~  r  T  *  “  7i!An 

3z  r  3r  3r  r  38 


(36) 


where  the  potentials  satisfy  the  homogeneous  scalar  wave  equations 
(25a).  (25b)  and  (25c) 


If  a  solution  where  there  is  no  8-dependance  is  examined,  it  is 
readily  seen  that  both  u^  and  u^  are  dependent  upon  both  ♦  and  A, 
while  Ug  is  dependent  solely  on  ¥.  Previous  models  [3)[l]  had  separate 
uncoupled  solutions  for  these  two  cases.  The  two-dimensional 
realization  based  on  ♦  and  A  could  support  compressional- ,  shear  in  r 
and  z  plane-  and  surface-waves  such  as  Rayleigh  ( vacuum- solid 
interface)  Scholte  ( liquid- solid  interface)  and  Stoneley ( solid- solid 
interface).  The  one-dimensional  realization  was  used  to  study 
transversal  shear  waves  (displacement  in  8-direction)  and  chiefly  Love 
waves.  These  are  transversal  shear  waves  confined  to  a  waveguide 
limited  by  two  interfaces. 


The  attenuation  in  the  medium  can  be  accounted  for  by  allowing  the 
Lame'  constants  and  thus  also  the  wavenumbers  to  be  complex.  This 
follows  from  including  viscous  losses  by  adding  a  first  order 
time-derivative  term  to  (1).  One  must  be  aware  when  including  losses 
in  this  manner  and  performing  integration  over  frequency  to  obtain  a 
time-dependent  solution,  the  answer  may  come  out  non-causal.  A 
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discussion  of  this  problem  is  found  in  [17],  This  point  will,  however, 
not  be  given  more  attention  here. 

In  order  not  to  confuse,  the  index  n  referring  to  the  layer  number 
will  be  implied  in  the  rest  of  this  section. 


A  set  of  solutions  of  equations  (25)  of  the  separable  form 


♦  ( r . 8 . 2 )  =  R(rl  0(81  Z(z)  (37) 

is  sought.  This  yields  the  seperated  ordinary  differential  equations 


0 '  ♦  m  0 


(38) 


2  2 

Z  *  (  h  -  k)  Z 


l  39  ) 


2  2  2  2 
r  R  .  rR  .  ( k  r  -  m  )R 


(401 


with  solutions 


0J (m0 )  =  COS (m8 ) 


(41a) 


(m8 I  =  sin (m8 ) 


(41b) 


Z  (  k  2  )  =  e 


-a(  k  I  z 


Zjlkz)  =  e 


al  k  )  z 


(42a) 

(42b) 


R  ,  (krl  =  J  (kr) 
m  t  m 


(  43a  ) 


**  -,(kr)  =  N  <kr)  l‘3b) 

mz  m 

respectively.  m  is  taken  to  be  0,1,2,...  as  the  boundary  conditions 

are  periodic  in  8 .  3  ()  and  N  (!  are  the  Bessel-  and  Neumann  functions 
m  m 

of  order  m  respectively.  The  Neumann  function  so]ution  (43b)  will  be 
discarded  as  it  has  a  singular  point  at  the  origin.  This  will  be  given 
further  d:  jssion  under  numerical  integration. 

The  solution  (37)  becomes  an  angular  expansion  of  the  potential 
functions 
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♦  ( r .  0 ,  z ) 


r  ♦m(r,z) 
m  =  0 


cos Im0  ) 
s in ( m0 ) 


(44  ) 


A  (  r  .  8  ,  z  ) 


l  A  Ir.z) 
m=0 


cos (m0 ) 
sin (m0 ) 


(45) 


Y(r . 0 ,z) 


where 


l  fm(r.z) 

i=0 


sin !ffl8  ) 
-cos (m8 I 


(46) 


4m( r , z ) 

=  J[a™(k)  e 

0 

za(  k  I 

♦ 

a"(k)  eza,k) 

Jk  J  (rk) 
m 

dk 

(47  ) 

Am(r.zl 

=  f[b™(IO  e 

0 

20<k) 

b"(k)  e2@(kl 

]  3  (rk) 
tn 

dx 

(48) 

Tm(r. zl 

=  JCc'J'lkl  e 

0 

zp(k) 

m  z0(  k! 

c2(k)  e 

]k  3  (rk) 
m 

dk 

(49) 

(k2  -  h2,W2 

.  k2 

>  Re { h  2 ) 

a  (  s  )  = 

3(h2  -  kV/2 

.  k2 

4  Re(h2} 

(50) 

(k2  -  g 2 , 

.  k2 

>  Re(g2) 

|3(sl  = 

j(g2  -  *V/2 

.  k2 

4  Re(g2) 

(51  ) 

The  displacement  components  u  ,u  and  u„  are  expanded  according  to 

z  r  o 


u  (  r  ,  8  ,  z ) 
r 


ug ( r , 8 . z I 


m 

l  um(r,z) 
m  =  0  r 


cos (m8 ) 
s in ( m6  9 


t 


C  Uplr.z) 
m=Q 


sin (m0 ) 
-cos (m8 I 


(521 


153) 
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u 

z 


(r.B.z) 


l 

m  =  0 


z 


cos (m8  ) 
s in ( m8 ) 


(54  ) 


If  we  derive  the  quantities  u  and  u  By  a  direct  application  of  the 

r  o 

operator  defined  by  (341,(35)  on  (44), (45  )  and  (  47  1,1  48  ),  we  would,  as 
we  take  into  account 


d_ 

dr 


J  (sr) 
m 


s J  -  ( sr ) 

nu  1 


m 

r 


J  ( sr ! 
m 


(55) 


obtain  results  with  different  orders  of  the  Bessel  functions.  It  is. 
from  a  numerical  point  of  view,  desirable  to  have  just  one  Bessel 
order  in  each  equation. 


If  we  form  the  linear  combination  of  (34)  and  135)  with  (52)  and  (53) 
as  sum  and  difference,  we  get 


t  <um 
m=  0  r 

cos (m8 ) 
s  in  ( m0  ) 

m  sin(m8)  , 
“  u8  -cos(m8) 

= 

,  ^  . 

J-  1 

♦  -  ,  a.  .  i 

3_ 

IT  ♦  —  1 

(  2- 

3r 

r  30 

3r  *  r 

30 

dz 

dr 

l 

r 


3_ 

38 


)  A 


=  l  {( 
m=  0 


3_ 

3r 


m  |  ^m  sin(m8) 
-cos (m8 ) 


m  j  ^m  cos(m0) 

s  in  (mB  ) 
r 


3,3  m  ,  m  sin(m8) 

♦  ~  I  A  ,  .  , 

3z  ,  -cos (m8 I 

3r  r 


(56  1 


when  we  use  the  notation  defined  in  (44)  -  (49). 


We  observe  that  the  partial  differentiation  operation  ^lth  respect  to 
8  applyed  to  the  expressions  for  ♦,  T  and  A:  (  4  4  )  -  (  49  ).  leads  to  a 
constant  m  (and  of  course  alteration  of  the  sine-  and 
cosine- functions  I  ,  so  that  the  partial  differentiation  operator 
becomes  equivalent  to  the  operator  defined  by 
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«fc 


-  )J  (sr) 


sJ 


.  1  sr  ) 


157  ) 


Hence,  we  have  only  one  Bessel  order  for  the  positive-  and 
nega tive - linear  combinations  respectively. 


The  trigonometric  expansion  given  by  (56)  will  of  course  be  the  same 
for  all  layers,  differing  only  by  the  expansion  coefficients.  It  is 
therefore  sufficient  to  meet  the  boundary  conditions  in  terms  of  the 
expansion  coefficients 


u  (  r  ,  z  )  *,  u.lr.zl 
r  o 


,  3_  -  m  , . m  -  3_ 

(  3r  *  r  ’♦  ‘  1  3r 


E  ,  y"1 
r 


o  z 


dr 


-  I  (58) 

r 


which,  when  ((17)-(i91)  is  inserted  becomes 


m,  ,  m  m,  ,  -za(s)  -  m,  ,  za ( s ) 

u(r,z)»u(r.z)  i  J [ -  a  ( s )  s  e  ♦  a  „  (  s  )  s  e 

r  -  0  o’  2 

m,  .  „ ,  ,  - zB ( s  I  -  m.  ,  . ,  ,  zp ( s  ) 

*  b( ( s  )  B ( s  )  e  »  b? ( s )  B I s )  e 


m  - zB ( s )  m  zB l s )  ,  ,  ,  ,  .. 

♦  c,(s)  s  e  ♦  c„(.)  s  e  ]  s  J  .(rs)  ds 

i  2  m  ♦  i 


(59) 


The  coresponding  development  of  u^  comes  considerably  easier  as  we 
observe  that  the  part  of  the  differential  operator  (36)  which  operates 
on  A  is  d2  -  (3/dz|2.  which  by  (25C)  becomes 


(60) 


so  that  we  can  write 
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m,  r,  m,  ,  ,,  -za(s)  ,  ,  za(sl 

u^lr.z)  =  Jf-a^s)  a(s)  e  *  a  2  ( s )  ct(s)  e 


m ,  ,  - zg( s )  m,  ,  zg  ( s ) ,  ,  ,  . 

*b  (s)  s  e  *  bis)  s  e  ]s  J  (rs)  ds 

1  2  m 


(61  ) 


The  following  boundary  conditions  at  the  horizontal  interfaces  involve 
the  stress  components  T  ,T^  ar|d  T0z'  these  are  expanded  like 
111),  115)  and  (16)  respectively,  use  of  Hooke's  law  leads  to  the 
following  expressions  for  the  expansion  coefficients 


T  (r.z) 
zz 


A V 2 (r.z)  *  2|i  7"  um(r,z) 
dz  z 


rr  m,  ,  ,  „  2  2,  -za(k)  ,,,2  2,  m  za(k) 

=  p/fa  (k)  12k  -g  )  e  *  (2k  -  g  )  a  2 ( k )  e 

Q 

-dm(k)  2  k  01k!  e‘Z@lkl  *  o'T  ( K )  2  k  0  f  k )  e2@lk))  k  1  (rk)  dk  (62) 
i  2  m 


T  (r.z)  .  T.  (r.z) 
rz  -  8z 


,  3  ,  m  ,  m, 

u(  ^  (ur(r,zl  ;  u9(r.zl 


.  o  m.  m, 

Cr~  ♦  )  w  {  r  .  2  )  ) 

dr  r  z 


gf [ ♦ a™( k ) 2ka ( k ) e  za{k)  ♦  a™( k ) 2ka( k ) e2a^ k ^ 
0 


‘b^lk,  (2k^-g^)eZ^kl  i  b* ( k ) ( 2k2- g2 ) e20 lk) 


-c’^lklkfllkls  ♦  c^fkJkfHkJe^^^Jk  J  ( k  r  )  dk 

1  2  m*  I 


(  B  3  ) 


It  snould  be  noted  that  the  use  of  the  linear  combinations  (59)  and 

(63)  in  the  boundary  conditions  rather  than  the  components  has  the 

effect  that  all  coefficients  relating  to  the  unknown  arbitrary 

functions  (a  ,a,,b  ,b  ,  c  and  c  )  will  be  independent  of  the  Fourier 
12  12  1  2 

order  m.  This  is  obviously  very  important  for  the  efficiency  of  the 


numerical  solution. 
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,Wi.  i„.r.  ..i.  <»•  *'•■*  *'  ■“*  “  pr“"'' 7 


Tm  tr,z> 

2Z 


A724m(r.zl 


-Xh2  Iu>le-Za(kl-a*tl<>e2a(l<,l  X 
0  ' 


I6M 


A  simple  compressional  source  placed  at  the  point  will 

produce  the  following  field  in  an  infinite  homogeneous  medium,  [<*] 


♦  (r.e .2) 

1 


s,  “  "  -a(k) |z-z  | 

- — k  Z  e  cosmt  8-8  I  f  [  J  (r  k)  - — - 

4 ttuj  „m  l  :  m  l  a(k) 

m  =  0  0 


] k J  (rk)  dk 
m 


(65) 


where  S  is  the  source  strength,  which  is  generally  complex  to  account 


for  the  actual  phase  of  the  source.  The  factor  e 


due  to  the 


expansion  of  the  exponential  function  in  a  Neumann  series  [9],  is 
1  ,  m  =  0 


(66) 


2  .  m  >  0 


If  more  than  one  source  is  present  within  a  layer,  the  contributions 
are  simply  added  to  yield 


♦  ( r  .  8  ,z) 


-  -  N  -a( k I |  z-z  | 

—  l  £  J[  I  S  cosm(  8  -  8  )J  ( r  k  )  3 - 1 

4  ttuj  „  m  .1  l  m  l 

m=0  0  i=1 


a  (  k ) 


]  k  J  Irk)  dl: 
m 

(67  ) 


where  N  is  the  number  of  sources.  The  potential  eir.S.z)  is  now 
expanded  like  (44)  and  the  coefficients  are  easily  obtained  as 

♦m( r . 8 . z ) 


e  “  N  cos(m8.)  -a(k)|z-z  | 

— 25  Jt  l  S  0 1  J  ( r  k)  1  k  1  (rk)  dk 

4nu)  1  .  i  sin(m8  m  i  a  k)  m 

0  i=  1  i 


(68) 


The  expansion  coefficients  for  the  Cartesian  displacements  components 
follow  from  ((52  ).  (  530  and  (54)) 
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z 


e 

m 

4  TTtJJ 


-N  COSlmQ)  a ( k ) I z  2  I 

f[  l  s  1  sign (z-z  )  J  (r  k)  ealK,|i  i1  ]kj  (rk)dk 

i  ,  „  ,  l  m  l  m 

0  1=1  sin(m8  ) 

l 


(69) 


u  (r  ,z)  *  ua (r.z) 
r  a 


-  N 

77..  Jts  z  s 


ittUJ 


cos (mB . ) 

1  J  (kr  )  S- 


-a  I k ) |  z-  z 


i=  1  1  sin (mB  ) 
i 


a(k) 


)  k  J  (rk)dk 
m*  1 


(701 


and  use  of  Hooke's  law  yield  the  following  expansion  coefficients  for 
the  stresses  involved  in  the  boundary  conditions 


(  r  ,  z  )  = 


turn 


-  N  cos (mB  ) 

I C  ( 2k  -gZ)  Z  %  1  J  > 

0  i  =  1  sinlmfl  ) 


-a( k ) | z-z  | 

S - 1. 

a(  k ) 


lk3  Irk)  dk  171) 
m 


_  m  .  .  _  m 

T  r.z  *  T  r.z  = 
rz  Bz 

pe  ”  N  cos (mB  )  .  .  .  . 

*  — —  J 1 2k  Z  S  1  sign(z-z  )  J  (kr  )e  2  i  ]kJ  (kr)dk 

-  iiriu  ;  .  l  .  ,  .  .  lmi  m*  1 

0  i=l  sinlmfl  I 

l 

(72) 


In  the  case  of  a  fluid  layer  eq .  (71)  must  be  modified  to 


_ro  .  . 

T  (r.z) 
zz 


Ah  e  ••  N  cos(m8  )  -a(k)|z-z 

— - m  Jr  r  S  1  J  (kr  I  s - 4 

Q  1=1  1  sinlmfl  1  m  1 

l 


a  ( k ) 


]  k  J  I kr )  dk 
m 


(73) 


Here  only  simple  compres s iona 1  sources  have  been  considered,  but  shear 
wave  sources,  involving  the  potential  A  or  T.  can  be  treated  in 
exactly  the  same  way.  leading  to  integral  representations  similar  to 
eqs  (69)  -  (72)  for  the  field  parameters. 


olution  technique 


2.5.1  Numerical  integration  of  the  Bessel  transform 

The  integrals  of  ((471,(48)  and  (49)1  may  be  written  in  a  simplified 
manner  as 


J  G(k)  J  (kr)  dk  (74) 

0 

The  Bessel  function  may  be  written  as  the  sum 

2  J  (zl  =  H* '  1  <z]  .  H( 21 (z)  (75) 

n  n  n 

where  H(,,()  and  hI2I(I  represent  the  Hankel  functions  of  first  and 
n  n 

second  kind,  order  n  respectively.  For  the  time-dependent  solution 
choosen  here  (exp(jivt)l.  the  Hankel  functions  of  first  and  second  kind 
represent  in-  and  outgoing  waves  respectively.  As  there  are  no 
backscatter  elements  included  in  the  model,  the  Bessel  function  may  be 
approximated  by  the  Hankel  function  of  the  second  kind.  This 
approximation  is  of  course  not  valid  in  the  region  between  source(s) 
and  the  center  axis.  Otherwise  the  validity  of  this  approximation  is 
discussed  in  [5],  The  Hankel  function  is  substituted  for  its  large 
argument  approximation  as  given  by 


h(2,(z,  =  1 

V  IT  Z 

) 1/2  e'3'2  -  f 

'I' 

1  76  ) 

We  will  evaluate 

the  integral  at 

discrete 

steps  according  to 

k  =  k  ♦  r>Ak 

Q 

n  =  0,1,2, 

.  .  .  . N-  1 

(  77a  ! 

r  =  rQ  ♦  mAr 

m  =  0.1.2, 

.  .  .  , N- 1 

(77b) 

AkAr  =  2  w  /  N 

(77c) 

Insertion  o'  (76)  and  (77)  into  174),  and  substituting  order  of 
integration  and  summation  yields 
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2.5.2  The  alobal  matrix  method 


Solution  of  the  problem  at  hana  implies  finding  the  values  of  the 
potentials  a.b  and  c  that  satisfy  the  pertinent  source  and  boundary 
conditions . 

At  each  interface  u  and  T  must  be  continuous.  At  solid/solic. 
z  zz 

interfaces  we  must  in  addition  require  that  u  ,u„,T  and  T„  be 

r  8  rz  8z 

continuous.  At  solid/liquid  interfaces  the  shear  stresses  must  vanish. 
We  express  these  boundary  conditions  with  (31),  (35),  136)  and  the 

expressions  for  stresses  obtained  by  the  use  of  Hookes  law.  These 
boundary  conditions  may  readily  be  reformulated  in  terms  of  the 
angular  expansion  coefficients  which,  when  collected  in  a  column 
vector  become: 


_m,  . 

F  (  r  .  z  ) 


uzlr,z) 

m .  ,  m 

u  ( r  ,  z  I  *u„ ( r , z ! 
r  8 

m,  .  m . 

u  Ir ,  z)-u_ I r , z) 
r  o 


Tzz  <r'zl 
_rn  . 

T  r,z)*T  (r.z) 
rz  8  z 

m  .  . 

T  ( r  ,  z ) -T  r.z) 
rz  8z 


(73) 


The  unknown  potentials  are  so  defined  that  the  depth  coordinate  z 
within  this  layer  is  0  at  the  upper  layer  interface.  The  boundary 
conditions  at  interface  n  which  is  the  lower  interface  of  layer  n  may 
thus  be  stated: 

Fm(r,z  )  ♦  Fm(r,z  )  -  F™  ( r .  0 1  -  Fm  (r,0)  =  C  ISO) 

nn  nn  n*t  n*l 

where  subscript  n  denotes  the  layer  number,  z  denotes  the  thickness 

n 

of  layer  n  and  the  terms  with  tilde  -  as  before  denotes  source  terms. 
Insertion  of  (611,  (59),  (82)  {(61)  in  the  fluid  case)  and  (63), 

interchanging  order  of  integration  and  summation,  reduces  (80)  to  a 
set  of  integrands  which  must  vanish.  Thus  we  are  left  with  a  set  of 
linear  equations  in  the  unknown  potential  functions  and  the  source 
contributions  which  becomes 
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A  ,  B  -A,  B  *  R  -  R  , 

n,l  n  n*  1  ,  u  n*  1  n*  1  .  u  n ,  1 

where  the  unknown  potential  functions  for  layer  n  are 


3m  I  k )  -- 
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For  a  solid  layer  the  matrix  A 


for  the  upper  interface  in 
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while  for  a 

liquid 

layer , 

the  matrix 

becomes 

A  (k)  = 
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The  matrix  for  the  lower  interface  of  a  layer  is 
multiplying  each  element  with  the  appropriate  exponential 


101  ) 


(82) 

layer  n  is 


(83) 


I  81  ) 

obtained  by 

as  given 


by  (611  -  ( 64 |  with  z= 


1  e : 
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A  (k)  = 
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The  source  contribution  vector  is 


c  N  cos(m0  | 

A  [  S  1  J  Ikr  ) 

liruj  i=1  1  sin(m8  )  m  1 
J. 


-a| z-2. 


r  - sqn ( 2-2  ) 

1  *  l 

-k/a 
k  /a 

M<2k2-g2J/ a 
2kpsgnlz-z  J 
-2kpsgn(z-z  j 


where  z  =  0  and  z 

n 


for  upoer  and  lower  boundaries  respectively. 


finaly  the  local  sets  of  boundary  equations  are  mapped  into  a  glubal 
set  of  equations 
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The  global  matrix  method.  as  opposed  to  matnzant  methodsOJ  is 
presented  in  [9]  . 
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2.5.J  numerical  considers t i on s 

A  commonly  known  problem  with  the  matnzant  method  on  which  FF  D  models 
haye  been  baaad[3],  is  encountered  when  dealing  with  thick  layers  in 
the  evanescent  region  of  the  horizontal  wavenumber  spectrum.  The  wave 
amplitudes,  represented  by  depth  dependent  exponential  functions  with 
positive  real  arguments.  may  attain  large  values.  It  is  of  course 
obvious  from  physical  considerations .  that  the  energy  content  cannot 
grow  beyond  finite  values,  so  this  situation  must  be  considered  as  a 
numerical  artifact.  Ceveral  techniques  are  utilized  to  remedy  this 
situation. 

The  matrices  A  and  A  are  made  dimensionless  by  dividing  the 
n  ,  i  n  ,  u 

stress-  and  pressure-related  coefficients  by  u/‘s  and  by  the 
horizontal  wavenumber  respectively.  en  denotes  the  density  of  an 
intermediate  layer.  This  will  ensure  that  the  coefficients  are  within 
the  same  order  of  magnitude. 

Each  layer  is  described  in  a  separate  local  coordinate  system  with 
origin  at  its  upper  surface.  This  will  ensure  that  the  value  of  the 
depth  does  not  exceed  the  layer  thickness. 

the  order  of  potentials.  as  defined  by  (82)  ensures  that  the 
coe f f l c l en t s  .  which  attain  high  values  due  to  the  above  mentioned 
exponential  functions.  come  close  to  the  diagonal  of  the  global 
matrix. 


These  remedies  will.  together  with  standard  pivoting  by  columns, 
ensure  that  the  solution  of  (88)  by  means  of  Gaussian  elimination, 
will  be  u  iconditionaily  stable[7J. 

The  above  statement  must,  however,  be  slightly  modified  when  we  have  a 
source  in  a  thick  layer  with  evanescent  propagation  conditions. 
Numerical  instability  due  to  the  nonvanishing  pertinent  row  of  the 
nghthand  side  of  (88).  may  occur.  This  problem  is  easily  circumvented 
by  introducing  dummy  -  in te r face s  gust  above  and  below'  the  source(s). 
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Coherent  sources  on  Both  sides  of  a  thick  layer  may  also  represent  a 
numerical  problem,  but  as  this  situation  is  regarded  as  being  less 
important  in  most  physical  applications.  it  is  not  considered  a 
serious  limitation.  It  is  of  course  possible  to  solve  the  problem 
separately  for  both  these  sources  and  superpose  the  solutions. 


35 


; ,  6  Numerical  example:  a  point  source  in  -free  space 

The  example  of  a  point  source  in  free  space  is  a  well  suited  example 
t'or  verification  of  the  model  as  the  correct  result  is  readily 
calculated  as  201og(r).  Free  space  means  a  water  layer,  containing 
source  and  receiver,  surrounded  by  upper-  and  lowe r - ha  1 f s pa ce s 
consisting  of  water  with  identical  parameters. 

This  example  is  however,  rather  a  challenge  for  the  FFP  model,  as  the 
integrand  easily  becomes  undersampled  due  to  the  branchpoint  arising 
from  the  square  root  in  the  denominator  of  !  7 1 )  Ithe  reader  is  also 

referred  to  (501).  Her».  the  integrand  is  sampled  at  8196  points 

-  6  - 1  -  l 
ranging  from  628  xlO  m  to  0.4  65  m 

Figure  2.1  shows  a  comparison  between  model  outputs  with  the  source 
located  on  the  center  axis  and  also  with  the  source  displaced  100  m  in 
positive  x-direction  from  the  center  axis,  l  e  the  model  is  run1 in 
two-  and  three-dimensional  modes  respectively.  In  both  cases  the 
frequency  is  100  He.  the  receiver  is  located  40  m  above  the  source. 
The  three-dimensional  case  required  an  angular  expansion  order  of  110. 
Curves  a  and  b  show  the  transmission  loss  for  the  two-dimensional-  and 
three-dimensional-cases  respectively.  The  region  of  validity  for  the 
computed  solution  begins  at  a  greater  distance  from  the  source  in  the 
three-dimensional  case,  as  the  approximation  to  the  higher  order 
Hankel  functions  has  larger  remainders  than  what  is  the  case  for  the 
0-order  approximation.  The  computed  solution  is  of  course  not  valid 
between  the  source  and  center  axis  as  only  che  Hankel - func t ion 
corresponding  to  outward  propagating  waves  is  included.  Figure  2.2 
shows  the  computed  transmission  loss  in  all  directions  for  the 
three-dimensional  example.  The  line  markers  x=0  and  y=0  are  drawn. 


The  objective  for  solving  the  inverse  problem  is  to  determine  the 
parameters  which  have  significant  influence  on  sound  propagation  in 
the  water  layer(sl  in  a  shallow  water  environment.  Determining  the 
parameters  within  '  acceptable '  confidence  bounds,  will  enable  us  to 
model  their  effects  and  thereby  predict  sound  propagation  in  the  water 
layer  under  the  influence  of  tne  ocean  bottom. 


To  my  knowledge,  there  has  not  been  any  earlier  reports  of  attempts  to 
infer  ocean  bottom  parameters  based  on  the  methods  to  be  utilized 
here . 


There  is  no  commonly  accepted  definition  making  a  clear  and  concise 
distinction  between  deep  and  shallow  water.  The  term  'shallow  water" 
is  used  here  to  imply  an  acoustic  environment  and  an  acoustic 
wavelength  sufficiently  large  (with  respect  to  oepth)  that  sound 
propagation  is  suited  for  modelling  by  wave  theory. 

While  sound  propagation  in  deep  wate:  is  well  handled  by  ray  theory, 
multiple  reflections  in  the  shallow  water  situation  create  a  rather 
complicated  ray  picture.  If  we, as  here,  limit  ourselves  to  a 
horizontally  stratified  medium,  multiple  bounces  may  interfere 
constructively  or  destructively.  and  thus  show  a  resonator  effect. 
This  effect  is,  as  will  be  demonstrated.  readily  described  in  the 
frequency-horizontal  wavenumber  IF-K)  domain.  A  projection  of  the 
signal  space  occuring  in  a  horizontally  layered  medium  into  the  F-K 
domain  will  readily  show  modes  of  propagation,  their  phase-  and  group 
velocities,  cut  off  frequencies  and  also  important  features  of  the 
ocean  bottom. 

The  term  mode'  is  used  in  a  more  free  sense  than  in  most  c  the  other 
literature.  Here  it  is  taken  to  mean  the  loci  of  points  in  the  =-K 
domain  where  constuctive  interference  occurs.  not  limited  to  rhe 
discrete  region  of  the  K-spectrum.  The  term  'normal  modes'  (normal  = 
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orthogonal!  does  not  enter  into  this  discussion  as  we  do  not  define  an 
appropriate  inner  product  onto  the  observation  space.  Elements  of  such 
an  inner  product  space  would  be  values  as  a  function  of  depth. 

The  nomenclature  for  characterization  of  interface  waves  is  adopted 
from  [  18)  .  A  Rayleigh  wave  can  occur  on  the  interface  between  a  solid 
and  vacuum,  a  Scholte  wave  between  a  solid  and  a  liquid  and  a  Stoneley 
wave  on  the  interface  between  two  solids. 

A  suitable  observation  space  for  projection  into  the  F-K  space  is  a 
horizontal  array.  In  this  work,  data  from  a  seismic  research  vessel 
towing  a  super  long  air  gun  array  (SLAG!  and  a  S2-element  hydrophone 
array  is  used,  not  because  this  particular  configuration  is  ideal  for 
our  purpose,  but  because  it  is  available. 


* 


It  is  not  intended  that  this  chapter  shall  he  an  exhaustive 
develpoment  of  waveguide  theory.  The  purpose  is  just  to  establish  some 
fundamental  properties  which  will  aid  us  in  understanding  how  to 
interpret  an  F-K  diagram.  The  complete  theory  is,  of  course.  covered 
implictly  in  the  previous  section  on  the  forward  problem. 


The  essential  elements  of  a  water  layer  bounded  by  surface  and  ocean 
bottom  are  shown  in  fig  ire  3.1.  Refering  to  equation  1641.  the  depth 
dependent  part  of  the  separable  solution  consists  of  an  up-  and 
downward  propagating  component.  The  boundary  conditions  are  equivalent 
with  reflections  ocouring  at  the  boundaries.  Constructive  interference 
between  a  and  ai  will  occur  when 
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where  R  ( tp )  e ^  ^ J  are  plane  wave  reflection  coefficients  and  1  is  a 
characteristic  angle,  either  incident  or  grazing,  the  subscripts  s  and 
b  refer  to  surface  and  bottom  respectively.  (89)  is  equivalent  to 


R  R  =  1  ( 90a ) 

s  b 

and 

0  •  0  ■  Oaz  =  (  1  -  n )  2  it  .  n  =  O.I.2.3....  (3.9b.' 

s  b 

! 9 0 a  1  implies  a  lossless  medium  and  totally  reflecting  boundaries, 
corresponding  to  a  resonator  with  an  infinite  0- factor.  'his  is  of 
course  an  approximation.  We  will  not  dwell  upon  its  validity,  as  the 
purpose  for  its  application  is  to  establish  a  f enomenological 

framework  for  a  basic  understanding  of  waveguide  properties.  !30b)  is 
a  resonance  condition.  As  an  a pp roxima t ion .  we  will  assume  the 
reflection  coefficient  at  the  surface  to  be  -1  for  all  ip.  The 

reflection  coefficient,  as  well  as  the  impedance,  for  the  ocean  bottom 

represent  a  neccessary  and  sufficient  description  of  the  influence  of 

the  ocean  bottom  properties  on  sound  propagation  in  the  water  layer 
for  a  horizontally  stratified  medium.  The  effect  of  more  than  one 
ocean  bottom  layer  may  be  collected  in  the  ocean  bottom  reflection 
coefficient.  but  it  then  becomes  a  function  of  both  frequency  and 
horizontal  wavenumber  as  reflections  are  related  to  incidence  angle 
and  propagation  between  boundaries  is  related  to  vertical  wavenumber. 
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3.3  A  -test  case 

Figure  3.3  shows  an  F-K  diagram  of  a  test  case  with  parameters  as 
shown  in  figure  3.2.  The  source  is  7.5  m-  and  the  receiver  is  14  m 
below  the  ocean  surface.  The  diagram  is  the  module  of  the  preweighted 
integrand  given  by  (78).  Radiating  from  the  origin  are  five  line 
markers  which  in  descending  order  of  phase  velocity  («j/k)  indicate  the 
loci  of  points  constituting  constant: 

compre s sional  wave  phase  velocity  in  the  ocean  bottom 
-shear  wave  phase  velocity  m  the  ocean  bottom 
-Rayleigh  wave  velocity  at  the  ocean  bottom  interface 
-compressional  wave  phase  velocity  in  the  water  layer 
-Scholte  wave  velocity  at  the  water/sulid  interface 
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Figure  3.2  Environmental  parameters  for  single  bottom  layer  test  case 


On  the  ordinate  axis  ksO  .  n  the  order  of  increasing  frequency,  the 
0..  1..  2..  3.  and  }ust  barely  the  <►  .  mode  may  be  identified.  The  0. 
mode  originates  at  the  origin. 


Figure  3.4  shows  module  and  phase  of  both  reflection  coefficient  and 
ocean  bottom  impedance  for  the  same  test  casp.  The  ipft  column  from 
top  to  bottom  shows  module  of  reflection  coefficient  including  the 
evanescent  region,  module  of  reflection  coefficient  excluding  the 
evanescent  region  and  finally  the  phase  of  the  reflection  coefficient 
including  the  evanescent  region.  The  right  column  shows  module  and 
phase  of  ocean  bottom  impedance,  both  including  the  evanescent  region, 
For  abscissa  the  horizontal  phase  velocity  is  choosen  as  it  offers  an 


with  the  characteristic  velocities  of  the  environment. 


easy  comparison 

The  horizontal  phase  velocity  relates  to  various  alternative 
parameters  as 
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f  where  u  and  are  incident-  and  grazing  angles  respectively  and  p  is 

‘  i  9 

slowness. 

Figures  3.5  to  3,10  are  identical  to  figure  3.4.  except  for  the  broken 
lines  which  indicate  per  (surba  tion  of  environmental  parameters.  In 
order  of  appearance,  the  effects  of  the  following  parameter  variations 
are  shown: 

- compres sional  wave  phase  velocity  m  the  water  layer 
-compressional  wave  phase  velocity  in  the  ocean  bottom 
-shear  wave  phase  velocity  in  the  ocean  bottom 
-density  in  the  ocean  bottom 

- compre s s iunal  wave  attenuation  in  the  ocean  bottom 
-shear  wave  attenuation  m  the  ocean  bottom 


The  source  and  receiver  are  located  fairly  close  to  the  ocean  surface. 
As  we  are  dealing  with  a  homogeneous  water  layer,  there  is  no  inherent 
mechanism  for  trapping  acoustic  energy  in  a  channel  close  to  the 
surface,  hence  energy  produced  by  the  source  and  detected  by  our 
receiver  will  be  influenced  by  bottom  properties  and  constructive 
interference  as  indicated  by  (89). 

The  0.  and  the  higher  order  modes  will  be  commented  on  separately,  as 
their  behaviour  differs. 


3.3.1  The  0 ,  mode 


The  0.  mode  originates  at  the  origin,  it  has  no  pronounced  low 
frequency  cut-off.  At  the  origin,  where  w  and  s  approach  zero.  the 
vertical  wavenumber  ((50). (26))  aproaches  zero.  As  the  depth  is 
constant,  az  in  (90b)  may  be  taken  to  be  zero.  thus,  for  n=0.  the 
phase  of  the  bottom  reflection  coefficient  approaches  n. 
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As  is  well  known,  the  Rayleigh  wave  occurs  at  the  interface  between  a 
solid  and  vacuum.  In  the  solid,  it  is  evanescent  in  the  z-direction 
and  is  thus  a  surface  wave  with  horizontal  wave  velocity  below  the 
free  shear  wave  velocity  of  the  solid.  For  a  liquid/solid  interface, 
The  Rayleigh  wave  becomes  a  pseudo- Rayleigh  wave:  its  horizontal  phase 
velocity  becomes  complex.  For  a  vacuum/solid  interface.  a  Rayleigh 
wave  implies  vertical  displacement  and  vanishing  normal  stress  at  the 
solid/vacuum  interface:  vanishing  impedance.  The  curve  for  module  of 
ocean  bottom  impedance  in  figure  3.1  shows  a  dip  at  what  is  here  the 
pseudo- Rayleigh  wave  horizontal  phase  velocity:  1  865.  S  m/s. 

From  the  graph  of  the  ocean  bottom  reflection  coefficient  in  figure 
3.1,  it  can  be  seen  that  while  the  module  has  a  dip  at  the  Rayleigh 
velocity,  the  phase  equals  -u.  Solution  of  ( 9 D b )  for  n  =  0  requires  that 
the  vertical  wavenumber  is  equal  to  zero.  Returning  to  figure  3.3,  we 
can  see  that  as  the  0.  mode  developes  from  the  origin,  at  a  phase 
velocity  asymptotically  equal  to  the  Rayleigh  wave  velocity  (actually 
figure  3.3  does  not  show  this  acuratly,  but  a  detail  of  this  corner 
near  the  origin  would),  its  energy  content  increases  due  to  increasing 
module  of  reflection  coefficient.  For  increasing  frequency,  figure  3.3 
shows  the  phase  velocity  to  decrease  and  go  below  the  water  velocity 
and  approach  the  Scholte  wave  velocity.  From  figure  3.1:  as  the 
horizontal  phase  velocity  decreases  from  the  Rayleigh  value.  the 
module  increases  and  the  phase  increases  through  a  peak  and  returns  to 
-a  at  the  water  wave  velocity.  Again,  we  have  a  solution  where  the 
vertical  wavenumber  a  is  zero,  but  now  because  the  horizontal 
wavenumber  equals  the  wavenumber  of  a  free  plane  wave  in  water.  In 
further  decreasing  the  horizontal  phase  velocity,  we  come  into  the 
evanescent  region,  characterized  by  a  real  (pure  real  for  a  lossless 
medium)  vertical  wavenumber,  As  we  enter  into  the  evanescent  region, 
the  acoustic  coupling  between  surface  and  bottom  diminishes,  and  the 
0.  mode  fades  out.  Obviously,  if  properly  excited,  the  0.  mode  will 
transfer  energy  into  a  pure  Scholte  surface  wave.  The  Scholte  wave 
velocity  is  for  this  case  1399.2  m/s.  Figure  3.4  shows  a  strong  peak 
in  module  of  reclection  coefficient  at  the  Scholte  wave  velocity. 
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Now  having  given  an  outline  of  the  Basic  mechanisms  of  the  0.  mode, 
we  can  use  the  perturbated  reflection  coefficients  of  figures  3,5  to 
3.10  to  look  into  how  variation  of  the  environmental  parameters 
manifest  themselves. 

As  is  well  known,  the  Rayleigh  velocity  is  governed  by  properties  of 
the  solid,  while  the  Scholte  root  is  governed  by  water  properties  as 
well.  For  more  details,  the  reader  is  referred  to  appendix  1  and  also 
to  [18]. 

The  location  of  the  dip  in  module  of  reflection  coefficient  at  the 
Rayleigh  wave  velocity  is  highly  dependent  upon  the  shear  wave 
velocity.  figure  3.7,  and  barely  upon  compre s s lonal  wave  velocity, 
density  and  shear  wave  attenuation.  The  dip  value.  however,  is 
strongly  dependent  upon  the  shear  wave  velocity  and  attenuation. 


3,3.2  The  higher  order  modes 


The  higher  order  modes  ( n  =  1  ,  2  .  -  .  in  (90b))  have  the  same  general 
behaviour  and  will  therefor  be  commented  upon  collectively.  It  is 
found  advantageous  to  divide  the  higher  order  modes  into  separate 
regions  characterized  by  the  horizontal  phase  velocities  as  follows: 


-continuous  region: 
-intermediate  region: 
-discrete  region: 


<  - 

<  c 

c 

<  c 

s 


As  mentioned  earlier  in  this  chapter,  the  F-K  plot  in  figure  3.3  is 
separated  into  these  regions  by  the  line  markers. 
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The  continuous  region.  For  the  continuous  region,  we  note  that, 
contrary  to  the  0.  mode,  the  higher  order  modes  have  a  low  frequency 
cutoff.  This  is  clearly  seen  in  figure  3.3  where  the  modes  originate 
from  the  axis  ssO.  For  s=0,  the  vertical  wavenumber  a  becomes  w/c. 
From  figure  3.4  and  figures  3.5  to  3.10,  it  is  seen  that  for  this 
region  the  phase  of  the  reflection  coefficient  equals  zero  for  all 
parameter  variations.  Consequently,  (90b)  simplifies  to 

"  z  :  *  1 2n- 1  )  .  n=  1  , 2 . 3 .  (92) 

c  2 

As  the  horizontal  phase  velocity  is  greater  than  the  compre s s lonal 
(and  of  course  the  shear)  velocity  of  the  ocean  bottom,  energy  is 
leaked  into  the  bottom  as  both  compressional-  and  shear  waves,  and  we 
do  not  have  a  trapped  wave  in  the  water  layer.  Consequently,  the 
transmission  loss  is  so  high  in  this  region  that  it  is  usually 
considered  as  non  propagating.  It  should  be  noted,  however,  that  for 
this  test  case,  we  are  considering  a  hard  ocean  bottom  with 
significant  impedance  contrast  between  water  and  bottom.  Hence,  even 
for  the  continuous  region,  the  reflection  coefficient  is  appreciably 
high . 

While,  as  a  function  of  increasing  horizontal  wavenumber,  the 
horizontal  phase  velocity  decreases  from  infinity,  the  group  velocity 
(dui/ds)  increases  from  zero.  For  this  region,  the  increase  in  group 
velocity  is  caused  solely  by  variations  in  frequency  and  horizontal 
wavenumber  such  that  the  vertical  wavenumber  a  remains  constant. 

For  the  continuous  region,  the  module  of  reflection  coefficient  is 
dependent  on  the  impedance  constrast  between  water  ano  solid,  hence 
the  phase  velocities  and  density.  As  the  phase  remains  unchanged, 
these  parameters  do  not  influence  the  loci  of  maxima  in  the  F-K 
domain , 

Horizontal  phase  velocity  equal  to  ocean  bottom  compressional 
velocity . 

The  highest  velocity  line  marker  in  figure  3.3,  indicating  a 
horizontal  phase  velocity  of  40QQ  mis  equal  to  the  compressional  wave 
velocity  of  the  ocean  bottom,  passes  through  peaks:  one  for  each  of 
the  higher  order  modes.  These  peaks  are  clearly  seen  as  a  peak  in 
module  of  reflection  coefficient  in  figure  3.4.  Figure  3.6  shows  how 
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the  peak  “follows"  ocean  bottom  comp  re ss tonal  wave  velocity.  Figure 
3.7  shows  how  its  shape  is  influenced  by  ocean  bottom  shear  velocity. 
It  should  also  be  noted  that  the  peak  is  slightly  influenced  by  ocean 
bottom  compres s ional  wave  attenuation,  unaffected  by  ocean  bottom 
shear  wave  attenation  and  also  that  it  is  approximatly  at  this  point 
the  phase  of  the  reflection  coefficient  begins  to  differ  from  zero. 

Intermediate  region. 

In  the  intermediate  region,  free  propagating  shear  waves  may  be 
coupled  into  the  ocean  bottom,  while  compressional  waves  undergo  total 
reflection.  The  group  velocity  increases  significantly  in  this  region. 

Horizontal  phase  velocity  equal  to  ocean  bottom  shear  velocity. 

When  the  horizontal  phase  velocity  equals  the  ocean  bottom  shear 
velocity,  the  bottom  becomes  totally  reflecting.  Hence  neither 
compressional-  nor  shear  waves  are  coupled  into  the  ocean  bottom  as 
free  propagating  waves.  According  to  [t],  the  group  velocity  equals 
the  phase  velocity;  the  ocean  bottom  shear  wave  velocity. 

The  discrete  region. 

The  discrete  region,  so  termed  because  the  horizontal  wavenumber 
spectrum  for  a  given  frequency  becomes  discrete  and  contains  the 
propagating  part  of  the  energy.  As  neither  shear-  nor  compressional 
wrves  are  transmitted  into  the  ocean  bottom  half  space,  all  energy  is 
trapped  in  the  water  layer:  we  have  a  wave  guide  effect.  In  this 
region,  the  phase  velocity  converges  from  above-  and  the  group 
velocity  converges  from  below  to  the  phase  velocity  of  the  water 
layer . 
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3.4  Description  of  data  collection  and  system  response 

The  measured  data  have  kindly  been  supplied  by  the  Norwegian  Petroleum 
Directorate.  They  were  gathered  by  the  seismic  essel  M/v  Malene 
Dstervold  operated  by  GECO.  The  source  was  an  airgun  array  and  the 
receiver  was  a  seismic  streamer,  both  towed  behind  the  vessel. 

The  data  were  collected  on  3t  July  1979,  north  of  Bear  Island.  75°?0’N 
0 

1  8  50  E.  This  area  is  characterized  by  its  hard  bottom,  thus  giving 
rise  to  multiple  reflections  between  ocean- sur f a ce  and  -bottom. 

The  data  are  identified  as  line  752Q-79  shots  nr.  9191  to  9296. 

The  streamer  is  towed  at  K  m  water  depth.  It  consists  of  52 
hydrophone  groups,  the  distance  between  their  centerpoints  are  50  m.  A 
group  is  formed  of  two  subgroups  each  of  length-  22  m  and  consisting  of 
32  equispaced  hydrophones.  The  spacing  between  subgroups  is  3  m.  The 
signals  from  each  hydrophone  within  a  group  are  summed  without 
weighting . 

The  hydrophone  group  separation.  50  m,  implies  a  spatial  angular 

sampling  rate  (wavenumber!  of  0.126  m  1  and  a  Nyquist  rate  of  D.D63 

- 1 
m 


When  taking  the  effects  of  the  ocean  surface  into  account  as  the 
Lloyd-mirror  effect,  the  modulus  of  a  horizontal  wavenumber  respons 
becomes 


Hst  =  ' 


—  (1  ♦  e35*  )  (1  ♦  e'32^  )  | 

.-37 


193) 


where 

T  =  s5  1941 
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V 


2.1/2  H 

s  )  d , 


(95) 


and  6  is  the  hydrofon  spacing  (22/31  m) .  1  is  the  distance  between  the 
first  hydrophones  in  a  subgroup  (25m),  N  is  the  number  of  hydrophones 
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within  a  subgroup  (32)  and  6^  is  the  nominal  streamer  depth  (Km). 

The  signals  from  each  group  are  sampled  successively  at  a  rate  of  250 
H 2 .  There  has  not  been  performed  any  filtering  in  the  low  frequency 
end  of  the  signal  spectrum,  the  high  frequency  end  has  been  filtered 
through  a  low-pass  filter  3  dB  down  at  64  Hz  decaying  72  dB/octave. 

The  distance  from  the  source  centerpoint  to  the  middle  of  the  first 
hydrofon  (offset)  is  180.75  m. 


The  source  has  been  towed  at  a  nominal  depth  of  7.5  m.  The  airgun 
array  was  fired  every  50  m.  The  frequency  spectrum  for  the  source  is 
shown  in  figure  3.12.  If  we  assume  the  source  geometry  depicted  in 
figure  3.11  and  include  the  Lloyd-mirror  effect,  we  obtain  the  modulus 
of  the  source  horizontal  wavenumber  spectrum  as 

H  =|(1*  e"')2r2)  t  e:  5iS  |  (96) 

30  1=1 


where 

r2  =  ((V  -  sV/2  d2 


(97) 


6  denotes  the  distance  from  first  to  l-th.  source  and  d  the  source 
l  2 

depth.  The  effects  of  source  spacing  transversal  to  centerline  of 
ship's  track  have  been  disregarded. 

Multiplication  of  (93)  and  (96)  gives  an  estimate  of  the  modulus  of 
the  horizontal  wavenumber  lespons  of  the  data  collection  system.  This, 
together  with  a  piecewise  linear  approximation  of  the  source  frequency 
spectrum,  is  shown  in  figure  3.13.  The  response  is  normalized  with 
respect  to  total  power.  Figure  3.13  a  shows  a  surface  plot  on  a  linear 
scale,  while  figure  3.13  b  shows  the  same  data  in  a  contour  plot  on  a 
logarithmic  scale  with  5  dB  between  contours.  Lloyd-mirro.  effect  of 
source  and  streamer  appear  as  a  null  radiating  from  the  origin  to 
{0.12,28}  while  th^  null  originating  in  {0,51}  is  caused  by  the 
streamer.  The  corresponding  null  produced  by  the  source  falls  outside 
the  figure  bounds.  The  source  sidelobes  are  seen  as  maxima  ar.d  minima 
at  constant  k-values.  The  streamer  produces  a  first  minimum  at  a 
horizontal  wavenumber  of  approximately  0.12  m~ 1  and  therefor  acts  as  a 
spatial  low  pass  filter  with  first  minimum  at  twice  the  Nyquist 
frequency . 
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3.5  Preprocessing  and  display  of  measured  data 

In  order  to  bring  the  measured  data  on  a  form  comparable  to  the  model 
output,  inverse  weighting  corresponding  to  the  post-FFT  weighting 
given  in  (78)  is  performed.  It  has  been  experienced  that  it  is  an 
advantage  to  Hanning  weight  the  data  in  the  X-direction. 

A  two-dimensional  FFT  is  performed.  The  FFT  lengths  are  64  points  in 
X-direction  and  1024  points  in  T-direction. 

The  displayed  data  is  limited  to  60  points  in  K-direction  and  250 
points  in  F-direction.  It  is  common  within  the  seismic  community  to 
display  the  K-direction  showing  both  positive-  and  negative  values  of 
horizontal  wavenumber,  but  as  we  have  good  signal  to  noise  ratio,  it 
is  fair  to  assume  that  energy  propagation  in  positive  direction  (from 
source  and  aft  along  array)  is  dominant.  This  agrees  well  with  the 
analyze!  data.  We  have  therefore  displayed  the  first  60  points  in 
K-direction  in  the  order  of  increasing  FFT-bin.  The  frequency 
direction  is  limited  to  250  points  as  we  do  not  have  significant 
energy  content  above  this  frequency. 

A  normalization  based  on  a  total  power  content  in  the  displayed  data 
to  be  equal  to  unity  is  performed.  The  amplitude  of  the  F-X  spectrum 
is  displayed  on  a  linear  scaled  surface  plot  and  a  logarithmic  contour 
plot  with  5  d8  between  contours  in  figure  3.14. 
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3.6  Modelling  of  the  data  collection  system 

Offset 

For  application  of  our  model. offset  is  defined  as  distance  from  axis 
r=0  to  first  receiver  position.  It  is  standard  input  and  follows  from 
equation  (78).  Offset  is  set  to  be  155.25  m,  see  figure  3.11. 

Source  array 

As  the  total  length  of  each  source  subarray  is  8.65  m  which  equals 
approximately  0.28  \  at  50  Hz.  we  will  approximate  each  subarray  as  a 
monopole.  The  source  array  is  readilly  modelled  by  application  of  the 
previously  derived  three-dimensional  capablities  of  our  numerical 
model.  The  source  array  is  aligned  with  the  axes  0=0,  z  =  7.5  m,  the 
center  subarray  on  the  axis  r=0. 

Receiver  array 

The  hydrophone  group  separation,  50  m,  sets  an  upper  limit  on 
observable  nonaliased  horizontal  wavenumber:  approximately  0.06  m  1  . 
Our  modelling  of  the  receiver  array  approximates  the  actual  array  by 
computing  a  4096  element  array  with  hydrophone  spacing  0.79  m.  This 
array  is  divided  into  64  subgroups,  each  subgroup  consists  of  64. 
hydrophones  with  spacing  between  hydrophone  groups  also  equal  to  0.79 
m.  Application  of  equations  (77)  with  hydrophone  separation  0.79  m 
gives  the  total  range  of  horizontal  wavenumber  to  be  8.855  m  1 .  It  is 
thus  sufficient  to  compute  a  horizontal  wavenumber  spectrum  ranging 
from  approximately  zero  to  0.12  m  1  and  set  the  region  ranging  from 
0.12  m  1  to  8.855  m  1  equal  to  zero 

As  previously  mentioned,  the  receiver  array  consists  of  52  hydrophone 
groups,  each  again  consisting  of  64  hydrophones.  The  approximations 
assume:  each  hydrophone  group  to  consist  of  64  equispaced  elements  as 

opposed  to  the  two  32-element  groups  spaced  3  m  apart  of  the  actual 
array,  that  the  separation  between  groups  is  0.79  m  as  opposed  to  3  m 
and  that  the  array  consists  of  64  as  opposed  to  52  hydrophone  groups. 
These  approximations  are  not  considered  to  be  of  significant 
consequence,  because,  as  demonstrated  in  figure  3.13.  the  role  of  this 
hydrophone  clustering  ls  to  create  a  spatial  low  pass  filter  with 
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first  minimum  at  approximately  0-12  m 

The  reasons  for  making  these  approximations  is  to  secure  that  the 
width  of  the  horizontal  wavenumber  bins  remain  unchanged  and  that  all 
energy  represented  by  the  original  horizontal  wavenumber  spectrum  is 
included  in  the  modelled  response. 

Th-  sequence  fur  modelling  of  receiver  array  is: 

-for  each  of  128  frequencies: 

-compute  64  point  horizontal  wavenumber  spectrum  from  0.-  to 
0.12  m'  1 

-preweight  according  to  (78) 

-include  zeroes  up  to  4096  point;  corresponding  to 
horizontal  wavenumber  9.355  m 
-perform  4096  point  FFT  over  preweighted  K-spectrum 
-post weight  according  t:  ■"9: 

-include  zeroes  in  frequency  bins  ’23  to  255 
-256  point  FFT  over  all  frequencies 
-  take  real  part 

-summation  over  64  hydrophone  groups,  each  with  64  hydrophones 
-inverse  postweighting  according  to  (78) 

-64  point  FFT  over  X-direction 
-256  point  FFT  over  T-direction 

Figure  3.14  shows  an  F-K  diagram  for  the  same  test  case  that  we  have 
investigated  previously,  except  for  that  the  source-  and  array 
configuration  is  modelled  as  described  here.  It  becomes  clear  that  the 
effects  of  source-  and  receiver  geometry  is  to  emphasize  the 
continuous  region  of  the  F-K  domain.  This  is  obviously  an  advantage 
for  the  se l smolog ica 1  purposes  for  which  the  system  is  designed,  as  it 
directs  energy  down  into  the  ocean  bottom. 

As  we  described  in  the  chapter  on  data  collection  system  response,  the 
effect  of  the  source  array  is  to  produce  a  downward  propagating  main 
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3.7  Geoacoustical  Background  information 


From  [19]  we  take  the  following  which  snows  the  prior 
geoacoustical  background  information  available  for  the  area  in 
our  shots  have  Been  recorded: 

water 


layer  number  1 : 

compre s s lona 1  wave  vel.  [m/sl 
density  I  kg/m3 1 
thickness  (ml 

layer  number  2: 

age 

sediment  i  rock  type 

compres s lona 1  wave  vel.  (m/s! 

shear  wave  vel .  (m/ s ) 

compre s s lona 1  wave  attenuation 

shear  wave  attenuation  (dB/M 

density  (kg/m3) 

thickness  (ml 

layer  number  3: 

age 

sediment  b  rock  type 

compres s lonal  wave  vel.  (m/s) 

shear  wave  •  el .  (m/ s ) 

compre s s lona 1  wave  attenuation 

shear  wave  attenuation  (dB/M 

density  (kg/m3) 

thickness  (m) 


1475-1185 

1000 

50-100 


gurassic-tnassLC 

shale-sandstone 

4200-4800 

2100-2400 

(d8/M  0.3-0.B 

1  .  0-2 . 8 
2450-2050 
800- 1000 


tnassic 
shale  sandstone 
5  100-5500 
2250-2750 

(d8/M  0.1-0. 5 

0.3-1  9 
2 .65-2.75 
800- 1200 


known 

which 
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layer  number  4: 

age 

sediment  &  rock  type 

compres s lonal  wave  vel.  (m/s) 

shear  wave  vel.  (m/s) 

compres sional  wave  attenuation  (dB/X) 

shear  wave  attenuation  (dB/X) 

density  ( kg/m3  ) 

thickness  lm) 


5800- 

2800- 

2800- 


1 


can  be  partly  covered  by  0.3-3  m  of  till,  glaciomanne  sediments  or 
Holocene  mud 
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.1 .  a  Identification-  and  modelling  of  events  in  the  measured  data 

F-K  diagrams  for  shots  numbered  1.  52,  100  are  shown  in  figures  3.15 
to  3.17  respectively.  As  the  sailed  distance  between  each  shot  is  50 
m,  the  distance  from  first  to  second  and  third  plot  is  2.5-  and  5  km 
respectively. 

Before  going  into  a  more  detailed  analysis,  we  would  like  to  note  the 
resemblance  between  the  three  figures,  which  of  course  suggests  that 
the  parameters  which  govern  propagation  in  the  water  layer  are 
basically  the  same  within  this  5  km  stretch.  As  a  consequence,  the 
environmental  parameters  we  may  identify  are  relatively  valid  for  at 
least  this  area.  It  is  in  all  figures  possible  to  identify  the  0.,  1., 
2.  and  3.  mode.  It  is  also  possible  to  identify  aliased  energy  in  the 
0.12  m  1  to  0.24  m  1  horizontal  wavenumber  region.  especially  in 
figure  3.17.  This  will  be  more  commented  on  and  exploited. 


3-8.1  First  approximation 

We  will  now.  having  established  a  framework  for  understanding  some  of 
‘he  features  of  the  F-K  plot,  proceed  to  attempt  to  interpret  the  F-K 
diagrams  for  the  three  shots  we  have  presented.  We  will  concentrate 
our  attention  on  the  first  shot,  shown  in  figure  3.15. 

We  will  start  out  by  determining  the  average  velocity  for 
compressional  waves  in  the  water  layer.  As  we  already  know  from  a 
previous  section,  the  phase  velocity  of  a  mode  will  approach  the  water 
velocity  assymptotically  from  above,  the  group  velocity  will  approach 
the  water  velocity  asymptotically  from  below,  and  hence  the  tangent  of 
a  mode  at  "high"  F-K  values  represents  the  water  velocity.  From  figure 
3.15,  a  trained  eye  can  observe  that  the  i.  mode  is  spatially  aliased 
and  thus  continues  from  approximately  (0,30).  The  contour  plot  of 
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figure  3.15  is  shown  in  figure  3.18  where  a  line  marker  for  K66  m/s 
passing  through  to.  12.28}  is  drawn,  also  for  the  spatially  aliased 
region.  The  first  mode  lies  above  the  line  marker:  its  phase  velocity 
is  greater  than  1166  m/s  and  its  group  velocity  (slope)  has  s  lower 
value  than  that  of  the  line  marker's  1166  m/s.  It  should  now.  if 
neccessary,  be  easier  to  identify  the  aliased  1.  mode  in  the  surface 
plot  of  figure  3.15. 

We  should  also  notice  that  the  1166  m/s  line  marker  crosses  the  0. 
mode  close  to  its  peak.  As  we  know,  this  peak  is  an  effect  of  both  the 
medium  and  the  source  geometry,  so  we  will  not  draw  any  firm 
conclusions  based  on  the  peak.  We  notice.  however,  that  the 
assymptotic  phase  velocities  are  above  the  line  marker  at  low 
frequency-wavenumbers  in  accordance  with  the  previously  described 
testcase.  The  0.  mode  will  be  looked  into  in  the  following. 

The  line  marker  in  figure  3.18  passes  through  (0.12.28}  in  the  F-K 
plot.  The  points  {0.12,27.5}  and  {0.12,28.5}  correspond  to  phase 
velocities  1440  m/s  and  1492  ml  s  respectively  (0.87.).  These  points 
are  in  the  order  of  reasonable  observation  accuracy.  Obviously,  this 
is  not  a  very  accurate  method  for  determining  ocean  sound  speed 
velocity  by  the  underwater  acoustician's  standards. 

We  now  proceed  to  estimate  the  cut-off  frequencies  for  the  1.,  2.  and 
3.  mode.  This  is  found  to  be  done  simplest  by  using  a  ten-point 
divider  along  the  frequency  axis.  With  point  0  at  the  origin,  point  1 
at  the  1.  mode,  point  3  at  the  2.  mode  and  so  on.  we  arrive  at  an 
estimated  cut-off  frequency  for  the  1.  mode  to  be  7  Hz.  Inserting  this 
value  and  the  phase  velocity  for  compressional  waves  in  the  water 
layer,  1466  m/s.  in  equation  (92).  yields  a  water  depth  of  52.3  m. 

The  surface  plot  in  figure  3.15  shows  two  very  distinct  peaks  in  the 
1.  and  2.  mode.  These  peaks  are  partly  a  result  of  the  source 
geometry,  but  they  are  much  sharper  than  the  lobes  of  the  source 
diagram,  see  figure  3.13.  We  take  these  peaks  to  be  at  the  phase 
velocity  of  the  compressional  wave  in  the  ocean  bottom  and  draw  a  line 
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through  them  as  shown  in  figure  3.19.  Again,  the  accuracy  is  limited. 
Out  if  we  accept  the  point  (0.1189,70)  (').  we  arrive  at  3700  m/s. 

Estimation  of  the  ocean  bottom  shear  wave  velocity  is  not  quite  so 
straightforward  as  the  other  quantities.  If  the  data  were  not 
" con t amina ted "  with  the  source  array  structure.  and  also  if  the 
hydrophone  spacing  were  closer  so  that  we  could  observe  a  higher 
horizontal  wavenumber  onaliased.  we  could  look  for  the  "point(s)' 
where  the  F-K  response  rises  up  and  becomes  discrete  and  also  perhaps 
for  the  pointlsl  where  the  group-  and  phase  velocities  are  equal.  As 
this  is  not  the  case  here,  we  shall  have  to  rely  more  heavily  on  the 
information  the  0.  mode  may  be  able  to  yield. 

An  estimate  of  the  Rayleigh-  and  Scholte  wave  velocities  for  the  ocean 
bottom  is  indicated  by  the  line  markers  in  figure  3.10.  It  must  be 
appreciated  that  the  determination  of  these  quantities  is  rather 
subjective  . 

The  Rayleigh-  and  Scholte  wave  velocities  for  som-  values  of  shear 
velocity,  computed  as  described  in  appendix  1,  is  shown  in  table  1. 
The  compres sional  wave  phase  velocities  for  water  and  bottom  are  those 
that  we  have  arrived  at:  1166  m/ s  and  3700  m/s  and  the  density  of  the 
ocean  bottom  is  assumed  to  be  2500  kg/m3 . 
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From  table  1.  we  see  that  for  shear  velocities  greater  than  the  water 
velocity,  it  is  the  p seudo - Ray le lgh  velocity  that  is  mostly  influenced 
by  variations  m  sound  speed.  All  the  computed  Scholte  wave  velocities 
lie  higher  than  1 2  0  m/s  as  estimated  in  figure  3.20.  We  will  for  the 
time  beeing  assume  the  shear  velocity  of  1915  m/s  which  matches  the 
estimated  pseudo - Ra yleigh  velocity,  1  7  8  0  m/s  figure  3.20.  The  line 
marker  for  1915  m/s  passing  through  C0.12.3G.57},  seems  to  be  somewnat 
too  low  for  the  group  velocity  to  equal  the  phase  velocity  in  the  1. 
mode  This  is  not  possible  to  determine  with  certamity.  but  it  would 
seem  that  a  somewhat  higher  shear  velocity:  in  the  Oder  of  1990  m/s. 
passing  through  {0.12,30}  m  figure  3.21,  would  give  a  better  fit  tu 
the  1.  mode.  This  implies  a  pseudo  Rayleigh  velocity  of  1814  m/s  (from 
table  1)  passing  through  {0.12.35.2}  which  is  acceptable  when  adopted 
to  shot  number  1  figure  3.21. 
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T he  geoacou s t lea  1  model  presented  in  chapter  3.7  suggests  higher 
compressional-  and  shear  velocities  in  the  first  bottom  layer: 
1200  -4800  m/s  and  2100-21U0  ml  s  respectively.  The  same  bottom  type, 
shale- sandstone ,  may  according  to  appendix  2,  however,  have 
compressional-  and  shear  velocities  in  the  range  2100-1800  m/s  and 
1200-2100  m/s  respectively.  The  estimated  values  should  thus  be  well 
within  geological  acceptable  Bounds. 

We  do  not  see  any  feasible  criteria  for  determination  of  ocean  bottom 
density  and  will,  based  on  the  geoacoustical  model  on  chapter  3.7, 
assume  it  to  be  2500  kg/m3  .  The  losses  will  not  be  considered  quite 
yet . 

Summarizing,  we  have  arrived  at  the  following  conclusions  for  a  first 
approximation  to  the  environmental  data  governing  shot  number  1.: 


compressional  wave  phase  velocity  in 

water  -. 

1466 

m/  s 

compressional  wave  phase  velocity  in 

bottom: 

37  00 

- 

shear  wave  phase  velocity  in  bottom: 

1990 

•• 

water  depth: 

52  .  3 

m 

density  of  bottom  (assumed): 

2.5 

kg/m' 

attenuation  of  comp  waves  in  bottom 

( assumed ) : 

2 

d8  IK 

attenuation  of  shear  waves  in  bottom 

( a  s  sumed ) : 

2  .  5 

- 

Line  markers  indicating  these  velocities  and  also  corresponding 
pseudo-Rayleigh-  and  Scholte  velocities  are  shown  in  figure  3.21 

The  modelling  of  this  1 .  approximation  to  shot  '  number  1  with 
simulation  of  source-  and  receiver  system  included  and  excluded  is 
shown  in  figures  3.22  and  3.23  respectively. 
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3.8.2  Second  approximation 


In  comparing  figure  3.22  with  figure  3.15,  the  most  striking 
discrepancy  is  to  be  found  in  the  continuous  region  wher  the  model 
predicts  a  much  higher  level  than  what  is  seen  in  the  measured  data. 

We  have  seen  from  figures  3.13  and  3.11,  that  the  main  influence  of 

the  source  and  receiver  geometry  is  to  emphasize  this  continuous 

region.  For  the  model  input  data,  we  have  appreciable  impedance 

contrast  between  water  and  bottom.  It  is  thus  in  no  way  unreasonable 
that  the  model  should  predict  this  response  in  the  continuous  region. 
We  must  therefore  look  into  what  effects  of  the  medium  we  have  not 
taken  into  account.  We  will  consider  the  following  hypothesis: 

-  the  source- streamer  geometry  is  not  according  to  specifications 

-  wrong  density  in  the  bottom 

-  wrong  attenuation  in  the  bottom 

-  sloping  bottom 

-  rough  bottom 

-  an  added  bottom  layer  with  high  attenuation 

An  erronuous  source- s t reamer  geometry  could  conceivably  be  such  that 
the  source's  mam  lobe  is  distorted  to  be  in  the  horizontal  wavenumber 
region  where  the  pronounced  spikes  occur,  approximately  0.01  m  1  to 
0.04  m  1 .  I  do  not  wish  to  speculate  on  the  likeliness  of  this 
hypothesis,  as  we  do  not  have  any  information  to  support  it. 

We  have  in  our  model  assumed  an  ocean  bottom  density  of  2500  kg/m" . 
Decreasing  the  reflection  coefficient  at  the  ocean  bottom  would  imply 
reducing  the  density,  see  figure  3.8.  The  absolute  lowest  acceptable 
value  is  1000  kg/m3,  otherwise  the  bottom  would  float  up  ('). 
Modelling  of  reduced  densities  have  shown  that  it  is  not  sufficient  to 
account  for  the  difference  between  model  and  measurement. 
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Figures  3.9  and  3.10  show  that  the  reflection  coefficient  in  the 
continuous  region  is  not  affected  by  variation  of  losses  in  the  ocean 
Bottom.  Therefore  the  hypothesis  of  wrong  attenuation  in  the  bottom 
cannot  be  accepted. 

Comparison  of  the  F-K  diagram  for  shot  numbers  1,  52  and  100,  figures 
3.15,  3.16  and  3.17,  show  that  they  grossly  resemble  each  other, 
consequently  a  steep  sloping  bottom  is  not  likely. 

For  the  frequency  regions  of  the  F-K  diagram  we  are  considering,  the 
wavelength  of  compressional  waves  in  water  is  in  the  order  of  50  m  and 
more.  As  we  have  arrived  at  a  water  depth  in  the  order  of  52  m.  it  is 
.tot  likely  that  there  can  be  bottom  roughness  on  a  scale  sufficient  to 
produce  the  reduction  in  continuous  region  response  we  are  looking 
for . 

Beside'  the  possibility  of  a  weakly  sloping  bottom,  we  are  left  with 
the  last  hypothesis:  an  added  bottom  layer  with  high  attenuation, 
inserted  "on  top  of”  the  existing  bottom.  We  should  at  once  be  able  to 
postulate  some  properties  this  layer  should  exhibit.  As  we  in  our  data 
“see"  the  bottom  we  have  arrived  at,  the  layer  should  be  transparent 
ir  the  sense  that  the  impedance  contrast  between  water  and  attenuating 
layer  is  almost  negligible.  Its  density  must  of  course  be  greater  than 
that  of  the  water. 

Conversion  of  energy  to  shear  waves  cannot  play  an  important  role  in 
the  continuous  region  as  the  particle  velocity  of  the  incident 
waterborne  wave  is  close  to  normal  to  the  boundary.  A  shear  wave  would 
have  to  be  downward  refracted  and  hence  have  its  'dominant  particle 
motion  parallell  to  the  boundary. 

The  lowest  line  marker  in  figure  3.23  is  set  at  1358  m/s  which  is  the 
computed  Scholte  wave  velocity  for  this  situation  (refere  to  table  1). 
The  0.  mode's  phase  velocity  is  above  this  value,  but  it  is  in  the 
process  of  reaching  it  asymptotically.  We  remember  from  figure  3.20, 
however,  that  we  have  estimated  the  asymptotic  value  for  the  0.  mode 
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to  be  at  most  1204  m/s.  For  our  situation  where  the  bottom  shear 
velocity  is  greater  than  the  compressional  wave  phase  velocity  in 
water,  it  is  this  water  velocity  that  has  significant  influence  on  the 
Scholte  wave  velocity.  We  will  choose  to  try  a  compressional  wave 
velocity  in  the  attenuation  layer  to  be  less  than  in  water.  for 
thereby  "pulling  down"  the  Scholte  wave  velocity  at  the  lower 
substrate  layer. 

The  geoacoustical  model  in  chapter  3.7  allows  for  a  0.3  -3m  thick 
layer  of  till,  gl3Ciomarme  sediments  or  Halocene  mud.  Appendix  2 
indicates  that  clay  exhibits  compressional-  and  shear  velocities  in 
the  region  100-2500  m/s  and  200-1000  m/s  respectively.  I  would  also 
like  to  add  that  I  have  participated  on  experiments  conducted  by  NDRE 
from  the  research  vessel  H  U  Sverdrup  in  the  same  area.  When  we 
recovered  equipment  from  the  ocean  bottom,  it  was  partly  covered  by  a 
blue  clayish  substance. 

Based  on  this  discussion,  we  will  introduce  a  thin  layer  with  high 
attenuation  and  acoustical  impedance  equal  to  water.  The  environmental 
parameters  for  the  2.  approximation  to  shot  number  1  will  be  as  shown 
in  figure  3.24. 
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The  losses  we  have  assumed  far  the  attenuating  layer  are  far  higher 
then  what  is  to  be  expected  in  known  bottom  types  [19],  but  they  have 
been  exaggeratet  in  order  to  see  if  attenuation  on  this  scale, 
whatever  physical  explanation  they  may  have,  could  account  for  the 
observed  low  response  in  the  continuous  region. 
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3.9  Discussion  of  fit  Between  measurement  ana  model 

F-K  diagrams  for  the  2.  approximation  with  environmental  data  as  shown 
in  figure  3.21,  including-  and  excluding  simulation  of  source - s treamer 
geometry  are  depicted  in  figures  3.25  and  3.26  respectively. 

Comparing  the  F-K  responses  with  and  without  attenuation  layer, 
figures  3.26  and  3.23,  we  can  see  that  the  attenuation  layer  has  had 
the  effect  of  reducing  the  continuous  F-K  region  response,  Put  as  can 
Be  seen  from  the  F-K  response  for  the  1.  approximation  when  we 
simulate  the  sou rce- receive r  geometry,  figure  3.25,  we  still  have  a 
far  too  strong  F-K  response  in  the  continuous  region.  Figures  3.27  to 
3.32  show  reflection  coefficients  and  Bottom  impedances  when  the 
parameters  of  the  attenuating  layer  are  perturbed.  We  notice  that 
shear  wave-  velocity.  attenuation,  and  density  do  not  affect  the 
module  of  reflection  coefficient  in  the  continuous  region. 
Consequently,  the  applied  shear  parameters  are  not  of  significant 
importance  in  the  continuous  region.  From  figure  3.32,  compared  with 
figures  3.27  and  3.30,  it  is  seen  that  on  a  relative  basis,  it  is 
layer  thickness  and  thereafter  compre s s lonal  wave  -velocity  and 
attenuation  that  have  influence  on  the  continuous  region  response. 

Again,  comparing  the  F-K  response  with  and  without  the  attenuating 
layer,  figures  3.26  and  3.23,  we  see  that  the  attenuating  layer  has 
recuced  the  spike  also  at  the  phase  velocity  equal  to  compres s ion  a  1 
wave  velocity  in  the  substrate  layer. 

Obviously,  we  have  not  reached  a  solution  with  a  perfect  fit  between 
measurement  and  model,  especially  on  a  relative  quantitative  scale.  I 
do  feel,  however,  that  there  is  some  success  on  a  qualitative  scale: 
the  water,  bottom  P-  and  S-  velocity  parameters  have  been  observed  and 
interpreted  in  the  measured  data  and  their  effects  in  the  sense  of 
loci  of  modes,  successfully  modelled,  Compr e s s lona 1  wave  velocities 
manifest  themselves  clearly.  The  shear  wave  velocities,  altough  they 
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have  a  more  obscure  manifestation,  are  collaborated  by  several 
indicators  (group  velocity  equal  to  phase  velocity  for  higher  order 
modes,  pseudo- Ray leigh  and  Scholte  wave  velocities  for  0.  mode). 

We  have  not  been  able  to  verify  that  the  assumptions  on  which  our 
numerical  model  is  based  upon,  are  fulfilled. 

The  question  of  uniqueness  of  solution  has  not  been  addressed  here. 
The  lobe  structure  in  the  measured  F-K  data  have  been  attributed  to 
the  side  lobe  structure  of  the  source  array.  Further  modelling,  not 
reported  here,  has  shown  that  such  a  lobe  structure  may  be  caused  by 
multiple  reflections  within  a  bottom  layer:  inter  bottom  layer 
resonance.  If  the  second  bottom  layer  has  a  greater  P-wave  velocity 
than  the  first.  P-waves  may  be  coupled  from  the  water  into  the  first 
layer  and  experience  total  reflection  at  the  boundary  between  1.  and 
2.  bottom  layer.  This  effect  can  manifest  itself  as  resonance  peaks  in 
tn-  region  where  the  horizontal  phase  velocity  is  between  the  P-wave 
velocities  of  1.  and  2.  bottom  layer.  This  hypothesis  has  not  been 
considered,  as  we  have  to  account  for  the  effect  of  source  sidelobes, 
which  are  in  reasonable  agreement  with  the  observed  lobes. 
Consequently,  a  combination  of  F-K-  and  time-  domain  methods  together 
with  a  monoploc  source,  may  prove  an  andvantage  in  solving  uniqueness 
problems . 

It  is  concluded  that  this  work  indicates  that  it  is  possible,  for 
situations  with  occurance  of  'many"  multiple  surface-bottom 
reflections,  to  infere  ocean  bottom  P-  and  S-  wave  velocities  from 
near  surface  measurements. 
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3.10  Suggestions  for  future  work 

For  future  work,  I  would  first  of  all  like  to  see  data  from  a  similar 
experiment,  but  with  a  different  source- streamer  geometry.  The  source 
should  be  a  monopole.  This  in  order  not  to  emphasize  the  low 
horizontal  wavenumber  region  and  also  to  avoid  sidelobe  contamination. 
This  would  readily  resolve  the  question  of  wether  the  lobes  observed 
in  the  measured  data  are  caused  by  source  sidelobe  structure  or  inter 
bottom  layer  resonance.  A  similar  subarray  of  the  type  used  in  the 
data  we  have  examined  here  should  suffice.  It  is  of  course  neccessary 
to  look  into  the  aspect  of  source  level  and  signal  to  noise  ratio.  For 
the  streamer  we  should  have  closer  hydrophone  spacing  so  as  to  be  able 
to  achieve  a  greater  unaliased  horizontal  wavenumber.  It  would  be  an 
advantage  to  repeat  the  experiment  in  the  same  general  area.  If 
another  area  were  to  be  chosen,  one  should  seek  a  flat  hard 
homogeneous  bottom  so  that  the  horizontal  stratification  assumption  of 
the  model  is  fullfilled  and  also  such  that  we  have  "many"  multiple 
reflections,  i  e  a  high  Q  resonator  which  consequently  shows  a 
pronounced  mode  structure. 

One  should  further  be  aware  of  the  water  depth,  such  that  the  number 
of  modes  is  within  practical  limits. 

One  should,  in  a  possible  continuation  of  this  work,  look  into  the 
feasibility  of  combining  time-  and  F-K-  domain  methods,  perhaps  in  an 
alternating  iterative  manner. 

Automatic  numerical  iteration  of  model  output  to  match  measurement  may 
be  undertaken  in  the  future,  but  initial  settings  of  the  model  must, 
in  my  opinion,  be  carried  out  prior  to  numerical  iteration. 


A  future  experiment  should  have  the  capabilities  of  ensuring  that  the 
assumptions  on  which  the  model  is  based,  are  fullfilled. 


A  numerical  model  of  wave  propagation  in  horizontally  stratified  media 
has  been  presented.  The  model  is  based  on  FTP- .  global  matrix-  and 
angular  expansion  techniques.  The  solution  is  exact,  except  for  a 
cylindrical  region  in  which  center  axis  and  sourcels)  are  included. 

Shot  data,  collected  by  a  seismic  survey  configuration  consisting  of  a 
source  array  and  a  streamer  is  presented.  Limitation  in  system 
response  imposed  by  source- streamer  geometry  is  discussed  and 
modelled . 

The  effects  of  environmental  parameters  on  events  observed  in  measured 
data  are  analyzed  and  modelled  in  the  F-K  domain.  The  fit  between 
model  and  measurement  compares  for  location  of  modes,  but  not 
completely  on  a  relative-quantitative  scale.  Uniquness  of  estimated 
environmental  parameters  is  discussed,  but  not  proved.  The  F-K  methods 
applied  here  should  be  supplemented  with  other  methods,  possibly  time 
domain  method s . 

A  future  experiment  with  monopole  source  and  closer  streamerhydro- 
phone  separation  should  be  carried  out.  It  shoul  preferably  take  place 
in  the  same  area,  or  in  a  similar  area  which  is  horizontally  statified 
with  multiple  reflections,  and  a  depth  so  that  the  number  of  modes  is 
within  practical  limits.  Core  samples  are  suggested  for  verification 
of  estimated  parameters 

The  methods  developed  here  indicate  that,  under  favourable  conditions, 
is  possible  to  infere  ocean  bottom  parameters  such  as  P-  and  S-  wave 
phase  velocities  from  near  surface  measurements. 
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Lunde  and  Trond  Jenserud.  Per  Riste  has  Been  of  great  help  in 
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involved,  is  also  highly  recognized  and  appreciated. 
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Appendix  1;  Rayleigh-  pseudo- Rayleigh  and  Scholte  waves 


From  D  Rauch  [18],  we  take  the  following  equation  for  determination  of 
Rayleigh-  and  Scholte  wave  velocities  at  a  boundary  between  a  solid 
and  vacuum  or  fluid  respectively. 


<X(X-,),/2,X-R,'/2 


(2X-1 ) 


IX-») 

(X-N) 


1/2 

1/2 


0  (A-  II 


where 


H  = 


X  = 


c 

s 

c 


and  the  subscripts  w,c,s  refer  to  water,  compres s lonal  and  shear 
properties  respectively  and  unsubscripted  c  is  the  the  horizontal  wave 
velocity.  For  the  case  of  vacuum  over  solid.  H.  and  consequently  the 
last  term  of  (A-1)  vanishes. 


For  both  cases  the  numerical  solution  is  found  by  stepping  X  through 
its  possible  range  and  seeking  the  minimum. 
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Jan  Inge  Faleide  [19]  has  compiled  the  following  tapel  for  some 
geoacoustical  properties  of  various  sediment  and  rock  type.  He  has 
collected  these  data  from  several  sources  to  which  he  referes 
collectively . 


P-VELOCITY 

S-VELOCITY 

DENSITY 

DENSITY 

(m/  s  ) 

1m/  s  1 

(g/cm3  ) 

1  aver ) 

AIR 

3 1 Q -  360 

0.0013 

OIL 

1250-1400 

0.60-0.90 

WATER 

1400-1550 

0.98-1 .05 

1  .025 

ICE 

3100-4200 

1600-2000 

0.88-1.07 

0.95 

CLAY 

1100-2500 

200-1000 

1 .50-2.60  ' 

2.20 

SILT 

1400-1800 

150-  450 

1 . 60-2 . 20 

1  .95 

SAND 

1000-2000 

100-  500 

1 . 60-2 . :o 

1  .90 

MORAINE 

1500-2700 

500-1300 

1 .50-2.00 

1  .  80 

SHALE 

2700-4600 

1500-2400 

2.00-3.20 

2  .40 

SANDSTONE 

2100-4500 

1200-2800 

2.10-2.80 

2.35 

CHALK 

2100-4200 

1 OOO-20OO 

1 .60-2.60 

2  .  00 

GYPSUM 

2000-3500 

1000-2000 

2.20-2.60 

2 .35 

ANHYDRITE 

3500-5500 

2000-3200 

2 .80-3 . 00 

2 .90 

SALT 

4200-5500 

2000-3200 

2.10-2.40 

2.15 

LIMESTONE 

3400-7000 

1800-3400 

2. 10-2.90 

2 . 55 

DOLOMITE 

3500-6900 

2000-3800 

2.40-2.90 

2  .75 

GNEISS 

3500-7500 

1700-3600 

2.40-3.00 

2  .75 

MARBLE 

3750-6950 

2000-3800 

2 . 60-2.90 

2  .75 

GRANITE 

4750-6000 

2400-3800 

2 . 50-2.90 

2  .65 

BASALT 

5500-6400 

2700-3400 

2.70-3.30 

3  .00 

GAB8R0 

6450-6700 

3400-3700 

2.70-3.30 

3  .  00 

ULTRABASIC 

7400-8600 

3700-4400 

3 . 00-3.40 

3  .  20 

ROCK 
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